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ABSTRACT 


The purpose of this thesis is to investigate the use of 
Conformal Subdomain Basis Functions (CSBF) in the Method of 
Moments (MM) solution of a thin wire scatterer. The effect of 
using CSBF on the computed current and the scattered field is 
investigated by formulating and coding the MM solution for a 
thin wire loop and comparing the computed results for various 
loop sizes to measured data and two other MM codes. 
Significant reduction in the number of segments (and computer 
memory requirements) are found for loops with circumferences 
of less than one to two wavelengths for plane wave incidence. 
From these results, it is concluded that the use of CSBF will 
Significantly reduce the number of segments required for the 


MM solution of a spiral antenna. 
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I. INTRODUCTION 


Numerical techniques for solving electromagnetic scattering problems using integral 
equations and the method of moments (MM) are well known. The physical problem, 
specified by Maxwell’s equations and boundary conditions, is reduced to an integro- 
differential equation over finite domains, and solved using a procedure referred to as the 
method of moments [Ref. 1]. The unknowns (usually currents) are represented 
by a series of basis functions with unknown expansion coefficients. The MM process 
generates a Set of linear equations that must be solved simultaneously. Until recently, 
these techniques have been limited to small (1 to 10 wavelength) geometries because of 
computer run time and memory constraints. With the development of faster computers 
with more memory, the MM technique has increasing application to larger geometries. 
However, computer memory and run time can still be inadequate to solve many 
important antenna and scattering problems. Numerically efficient solutions require less 
computer memory and/or less computer run time. Therefore, any increase in the 
efficiency of a MM solution is of great practical interest. 

The usual MM approach to modeling a thin wire of arbitrary shape is to 
specify a series of points, with piecewise linear segments between the points to 
approximate the wire. The current is represented by one or more basis functions, each 
with constant phase, over a piecewise linear segment. Typically, the size of the segments 


is set by how accurately the current or scattered field needs to be determined. For 


convergence of the current, segment lengths of 0.05 A to 0.1 X are generally required. 

A second factor that influences the segment size 1s the radius of curvature of the 
wire. Tightly curved wires require smaller segments to reproduce the wire shape 
accurately. When the radius of curvature is larger than a wavelength, the first case sets 
the segment size; when the radius of curvature is much less than a wavelength, the 
second case dictates the segment size (Figure 1). 

All generally available MM codes based on the method of subdomains use 
the first approach. A natural question arises: Does dividing the wire into curved segments 
that conform to its shape, with arclengths restricted only by the maximum current 
variation rule, yield a solution that converges with fewer subsections? If the answer is 
yes, 1s the improvement in convergence worth the greater complexity and coding effort? 


To resolve this issue, the following approach is taken: 





Figure 1. Left: Piecewise Linear Segments, Right: Curved Segments 


1. Formulate the solution for linearly and circularly polarized plane wave incidence 
for a geometrically simple shape such as a loop. 


2. Computer code the solution in FORTRAN. 
3. Validate the solution using other MM solutions and measured data. 


4. Study the convergence of the solution with respect to loop parameters and 
compare its performance to a method using linear subsections. 


5. Study the effect of changes in program structure on its computational efficiency. 
Chapter II discusses the derivation and the MM solution of the electric field integral 
equation for a thin wire loop. Chapter III discusses three MM FORTRAN programs used 
to determine the current on the loop. The entire domain solution (due to R. F. 
Harrington), which uses complex Fourier modes (HARLOOP) is considered to be the 
most accurate and therefore serves as a baseline for evaluating the other solutions. A 
second program that uses linear segments (LOOPSCAT) 1s compared to a third program 
that uses curved segments (CURVENEW). Chapter IV discusses the results obtained by 
the three methods and presents some guidelines in choosing an optimum solution method 


for a given antenna or scattering problem. 


II. THE THIN WIRE INTEGRAL EQUATION 


A. DERIVATION OF THE THIN-WIRE ELECTRIC FIELD INTEGRAL 

EQUATION 

In this chapter, the integral equation for the current on a thin wire will be 
developed. Time-harmonic field quantities are assumed throughout. Phasor quantities are 
used with the et dependency suppressed. 

Referring to the thin wire geometry of Figure 2, the origin 1s point QO, the location 
of a source point is given by the vector r' and an observation point by the vector r. The 
wire radius, a, 1S considered constant over the length of the wire. The vector | 1s 
everywhere parallel to the surface of the wire. Since the sum of the tangential 
components of the incident and scattered electric field must vanish at the surface of a 


perfect electric conductor, the boundary condition is, 


i x (E'+E‘) = 0 (2.1) 


If the radius of the wire is small compared to the wavelength of the excitation, the 
surface current density, J,, can be considered constant around the circumference of the 
wire and directed along its axis. The excitation field can be either an incident wave or 
an impressed voltage. The scattered field is the field due to the current on the conductor 
induced by the excitation field. 


The wave equation in terms of the vector potential A is given by 


V *A+B7A = -pJ, (2.2) 
where 8 =27/X. The solution to equation (2.2) is 


J oe 
ie - f+~—_—as’ = pfJ, g(r,r’)dS' ae 
N's s 
where the integration is over the primed (source) coordinates. The Green’s function, 


g(r,r‘) is defined as, 


“iB \r-r'| 
eS _ (2.4) 


4n|r-r’| 


The expression for the scattered electric field in terms of A is, 





Figure 2. Thin Wire Geometry 





ES SOAS Ae (2.5) 
WPe 


Applying the boundary condition of equation (2.1), 





i $s : Jj 
B= oe = JwA+ V(V-A on S. (2.6) 
me Me Tl Ae (V-A) 


Substitution of A in equation (2.3) into equation (2.6) gives, 


re = jop[J,g(r.r')ds'« — V-: 





on S. (29) 








p[J,g0r')dS' 


S 


Call the second term on the nght side of equation (2.7) V, and assume the medium to 


be homogeneous, 











V = V1 Ree vad aS ev jul Ven iec(rmm)ae le (2.8) 
s’ 5. 
The vector identity for the divergence of a scalar u times a vector v 1s, 
V- (uv) = Vu-ve+u(V-y) . (2.9) 
Applying this identity to equation (2.8) gives 
V=Vin[(Ve@r’))-J,ds/| . (2.10) 








5s 


It can be shown that Vg(r,r') = -V'g(r,r') [Ref. 2]. Using this in equation (2.10) 


and applying the identity of equation (2.9) again yields, 


ee pf (Vg (rr'))-J,ds’ 


s 


pf gcnr')W'-J,)d8!-nfV!> Gor!) J,)d5' 
Ss s! 








(2.11) 
Ay 





where V' is the del operator defined with respect to the primed coordinates. The second 
integral on the right side of equation (2.11) is equal to zero by the surface divergence 


theorem [Ref. 3]. Thus, V simplifies to 


Van [ (VT) Ve(r,r')d5! . (2.12) 


s! 


Substitution of equation (2.12) into equation (2.7) gives an integral equation for J,, 


B= = jou {J,g(rr!\ds! + [WT )Vg(rsr')d5' (2.13) 
s! 


s! 


which may be expressed more compactly as 





| See = [ jaws eter) + Led eetrr’) aS! (2.14) 
WE 
S 


f 


Equation (2.14) is a form of the Electric Field Integral Equation (EFIE). The unknown 


quantity to be solved for is J,. 


B. SOLUTION OF THE EFIE USING MM 
The method of moments (MM) technique can be used to solve for J, by expanding 


it into a series of basis functions, J,, 


N 
=) Ch (Za) 


where the C; are complex constants to be determined. Substitution of equation (2.15) into 


2.14 gives 


C. 


i 
Ss; 


4 


jops.g(rr)+L_(V'-J,) Vern’) | dS! . (2.16) 
WE 





eo, 
W 
jot 


To generate the required N equations to solve for the N unknowns, define a suitable 
weighting function W,, and take the inner product of W, with both sides of equation 


2.16. The inner product is defined such that it satisfies 


<w,y> = <v,H> 

<af+yv,w> = a<f,w>+y<v,w> (2.17) 
<y",w> >0 ify + 0 

<v",w> = 0 ifv =0 


[Ref. +]. Choose the following inner product: 


<W,y> = [words (2.18) 
§ 


where * signifies the complex conjugate. This leads to 


[W -E' dS = Liou, Le/| Jgnr')+L(v J) Vg(r,r')| dS'dS (2.19) 


5, 5, 5; 


Interchanging the order of summation and integration and applying the surface divergence 


theorem yields for the right hand side of equation (2.19) 


reff 


jon(W, -J,Je,r)--L(V-S,(V- W, grr’) |d5/dS (2.20) 
t=1 ScS, WE 





By making the following definitions, 


ali WeErexdS 
*t . (2.21) 
Zee i [ jop(W,-Jg(r,r')- (VI (V: W, )g(r,7')| dS'dS 
WE 
58 





2.20 and 2.21 can be wnitten in matrix form, 


(V] = (Z}[C] (2.22) 


where [V], [Z] and [C] are called the generalized voltage, impedance and current 
matrices, respectively. The unknown [C] may be solved by an appropriate matrix 


inversion algorithm. Symbolically, 


(Qe WAP a C23) 


The generalized current matrix elements are the weighting coefficients in the summation 
of equation (2.15). The current J, 1s computed from equation (2.15) and the scattered 
field is calculated using this current in the radiation integrals. It should be noted that [V], 
[Z], and [C] have units of volts, ohms, and amperes, but are not unique. In general, they 
depend on the choice of basis and weighting functions. However, the current will 
converge to the same numerical value as the number of basis functions are increased, 


provided the solutions are implemented correctly. 


C. SPECIALIZATION OF THE EFIE TO A CIRCULAR LOOP USING 

CONFORMAL SUBSECTIONS. 

The MM procedure will now be applied to a circular loop in the X-Y plane as 
illustrated in Figure 3. The loop is an ideal test geometry to study the characteristics of 
a MM solution using conformal subsections. It is a relatively simple geometry and other 
accurate solution methods are available to evaluate the performance of conformal 
subsections. The loop has a radius f, and is divided into N conformal segments. In this 


case a conformal segment is a circular arc. The arclength of the i" segment is, 


Al. =I! 


i i+] 


-1 = rAd, . (2.24) 


The basis functions, J,, of equation (2.15) are chosen to be overlapping triangles, 





Figure 3. Thin Wire Loop Geometry 
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HII) 











; eS Se 
2naAl, 
fee LUNE St TAN Il (2.25) 
[fis eee 
27a I & 


0; elsewhere 


Triangular basis functions are chosen because they are a more accurate representation of 
the current than a pulse basis function since the current is continuous everywhere along 
the wire, and they are relatively easy to deal with analytically. Balanis [Ref. 5] 
states that increasing the basis function complexity beyond triangles may not be 
warranted by the additional improvement in convergence rate. The triangular basis 
functions span two segments and overlap as shown in Figure 3. Therefore, the resultant 
current will be piecewise linear. The weighting (or testing) functions W, in equation 
(2.19) are chosen such that W,=J, (Galerkin’s method). Wang [Ref. 6] states 
that Galerkin’s method provides numerical results which are more accurate than other 
testing methods under similar computational constraints. Substitution of the above 
weighting and basis functions in the expression for Z, in equation (2.21) yields 


20 ff 


jopT (DT. )E-1) - a TU) T' ) \g(,r')aS'as 
S,S, 





(225) 


oT. oT 
Wiereeiua()) p= ee TO) c= 
al’ ol 


Using the following relations, 


ikl 


Guy 
mo) 
{I 


T= $-$' = cos(o -$’) 
rod; I’ = ro’ 

dS = 2nar,dp; dS’ = 2nar,do’ 
T.() = T(rod) = T,(2) 


be] 
i] 


TW) = —7',() (2.27) 
0 
n = [ B = w/pe 
As e FBIR| Le 
g(r,r’) = 4n|R\’ r-r 


equation (2.26) may be written as 





2p. O4-2 %i-2 _jB'R| 
Zip OE [ {| m@reoreoso-4')-— 7 6)7,(@)|£—doag! (2.28) 
an on. B?r° |R | 





From Figure 4 and the law of cosines, | R | is given by, 


|R|? = 2r,°[1 -cos(p-')] +a? (2°29) 





Figure 4. Geometry for Determining R. 


2 


and using the trigonometric identity, 
sn'() = = (1-cosa) (2.30) 


meSUIES iN, 


(225) 





By choosing the test (unprimed) points at the center of the wire and the source points on 
the surface of the wire, the singularities along the diagonal of [Z] at 6=¢' where r=r' 
are avoided. The technique used to calculate | R_| is discussed further in Chapters III 
and IV. 


The voltage elements, V,, given in equation (2.21) become 


To%4.2 
Vere [ T, (roo) Eid . (2.32) 


rod; 


The incident field, E', for the purpose of this study, will be a plane wave. Figure 5 
shows the direction of incidence of the plane wave in spherical polar angles 6=© and 
@= measured from the Z and X axes, respectively. E' can be @ or ¢ polarized. 


Referring to Figure 5, for 6 polarization, the component of E tangential to the loop, is 
E,®-$ = E,cos® sin(®-@) . (2.53) 


Similarly, for @ polarization, 





f 
x 
Figure 5. Plane Wave Incident on Circular Loop 
E.&-> = E,cos(®-4) . (2.34) 
The component of the phase vector, 8, parallel to r 1$ 
6 -# = sin® cos(®-¢) . (2.35) 
Equations (2.30) and (2.32) combine with 2.29 to give 
$4.2 
Vs, = ToEgcos® i T,()sin(® - p)e Mr OO ag (2.36) 
>} 
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A similar expression for a ¢ directed incident field 1s 


$4.2 
Ven = 1oBy | Ty) c0s(@- hye Presmerre ag, (2.37) 
$y 


The computer coding of the solution for the thin wire loop using the equations 


developed above is described in the next chapter. 
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Iii. COMPUTER CODES FOR THE THIN WIRE LOOP 


In this section, the FORTRAN program for a thin wire loop using curved segments 
is discussed. The results are presented and compared to similar solutions using straight 


subsections and Fourier modes. 


A. DESCRIPTION OF THE CODES 

Table 1 summarizes the organization of the three programs. Computer listings are 
given in Appendix B and equations from Chapter II will be referenced with a "2." 
preceding the equation number. The FORTRAN source code for the conformal 
subsections 1s named CURVENEW, and the codes for the straight subsections and the 
Fourier mode solution are named LOOPSCAT and HARLOOP respectively. 

CURVENEW, LOOPSCAT and HARLOOP are functionally similar. Each 
calculates the loop geometry based on the segment size, loop radius, and wire radius, and 
each uses Gaussian quadrature for numerical integration. CURVENEW computes the 
impedance matrix, [Z], 1n subroutine ZCURVED from the loop geometry of Figure 5 
using equation (2.28). LOOPSCAT uses a similar formulation applied to straight 
segments in subroutine ZMATWW. HARLOOP computes [Z] using the equations 
developed in reference [7]. CURVENEW computes the excitation vector, [V], using 
equation (2.36) or (2.37) in subroutine CURVEW. LOOPSCAT uses a similar 


formulation for straight subsections in subroutine PLANEW. HARLOOP computes [V] 
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using the equations in reference [7] in subroutine PLANEW. In CURVENEW, all 
integrals are evaluated numerically and symmetry of the impedance elements is used to 
fill the [Z] matrix and reduce the number of numerical calculations. Two formulations 
were investigated with LOOPSCAT: One using a delta function approximation to 
evaluate one of the double integrals in equation (2.21) and the other using Gaussian 
quadrature for both integrations. CURVENEW does all numerical integrations using 


Gaussian quadrature. Matrix symmetry 1s also used in LOOPSCAT to reduce the number 


TABLE 1. FUNCTIONAL SUMMARY OF PROGRAMS 


CURVENEW 
(Curved (Straight (Fourier 
Segments) Segments) modes) 


Lines 25-38 Lines 15-27 Lines 19-30 


Lines 45-7] Lines 53-76 Lines 44-54 


Program-- > LOOPSCAT HARLOOP 


Function 


Read Input 
Parameters 


Establish 
Loop 
Geometry 


Calculate 


[Z} 


Calculate 
[V] 


Solve 
System 


[C]={Zy'[V] 


Calculate E> 


Subroutine 
ZCURVED 


Subroutine 
CURVEW 


Subroutines 
DECOMP and 
SOLVE 


Lines 140-180 


Subroutine 
ZMATWW 


Subroutine 
PLANEW 


Subroutines 
DECOMP and 
SOLVE 


Lines 160-197 


Subroutine 
ZMATWW 


Subroutine 
PLANEW 


Subroutines 
DECOMP and 
SOLVE 


Lines 93-137 





of calculations for [Z]. Computation of [C] 1s performed in subroutines DECOMP and 
SOLVE [Ref. 3], which solve the system of equations using Gaussian elimination. 
Subroutines DECOMP and SOLVE are common to all three programs. The subroutines 


ZCURVED and CURVEW are discussed in more detail in the next section. 


1. Loop Geometry 
To generate the loop geometry an initial estimate of the desired circular arc 


length, Al, is provided. This 1s used to estimate an angular increment, Ad, 


Ad = —. (Sa) 


From this estimate, the number of generating points 1s calculated by, 


Nee el (3.2) 
Ad 
and Ad is recalculated using, 
27 
Nao (323) 
’ in 


to ensure that A¢ is such that N segments fill exactly 2a radians. The loop points P, and 
P, are coincident with P,, and P,,, so that the current is continuous around the loop. 
2. Subroutines ZCURVED and ZMATWW 
Subroutines ZCURVED and ZMATWW take advantage of the symmetry that 


exists on the loop. For all basis functions, the self impedance terms are equal. The 
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“remainder of the matrix is filled with impedance values that repeat in a cyclic manner. 


Mathematically, 
Z, = 2 = 233 = Ly 
Zi, = £53 = Ly, == Lyin 
TE 235 = £y-2N (3.4) 
Cari = ZyN 


The elements along any diagonal of [Z] are equal and the lower off-diagonal elements 
are the mirror image of the upper diagonals. Thus [Z] is a symmetrical Toeplitz matrix. 
Therefore, computation of the first row of [Z] provides enough information to fill the 
entire matrix. 

Because of the Green’s function in the integrand for the impedance elements, the 
numerical treatment of the self terin is very important. To optimize the convergence rate 
and accuracy of CURVENEW and LOOPSCAT, several different approaches are used 
to evaluate | R| near the singularity point where r=r'. In the first method, the 
observation point is chosen along the axis of the wire and the source point along the 
surface for all i,k, giving | R | =a at 6=¢' (equation (2.31)). For the second method, 
both the observation point and the source point are chosen along the axis of the wire 
except on the segment i=k where the value of ¢ at the midpoint is chosen on the axis 
of the wire, with r’ on the surface of the wire. Finally, both the observation point and 
the source point are chosen along the axis of the wire, except on the segment 1=k, where 


r is chosen along the axis, and r' is chosen on the surface. Choosing the source point and 


Jee, 


observation point as in the first case gives the most accurate results, but only slightly 
more accurate than the third case. Case two is accurate for small segment sizes but is 
inaccurate for larger segment sizes. Case three was Selected for both CURVENEW and 
LOOPSCAT because it is only slightly less accurate than case one, and required fewer 
lines of code. 
ZCURVED calculates the impedance elements of the first row of [Z] by breaking 
the integral in equation (2.28) into four parts. For example, in the first row of [Z], Z,,, 
is calculated by summing contributions from the following four regions of integration 
(Figure 6): 
1. A double integration along the positive slope of the T, over segment 1 and the 
positive slope of T; over segment 1. 


2. A double integration along the negative slope of the T, over segment 2, and the 
positive slope of T, over segment 1. 


3. A double integration along the positive slope of the T, over segment 1, and the 
negative slope of T, over segment 1+ 1. 


4. A double integration along the negative slope of the T, over segment 2, and the 
negative slope of T,; over segment 1+ 1. 
The integrations are computed in a similar manner for the derivatives of the basis 
functions, T,' and T;', over the same subsections. A similar procedure is used for the 
Straight subsections in subroutine ZMATWW to calculate the impedance elements. 
The numerical integrations are performed using Gaussian quadrature, with 
the number of points per subsection specified as an input parameter to program 


GAUSWGT, which computes the Legendre polynomial coefficients for a specified 
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number of integration points and writes them to file OUTGLEG. Gaussian quadrature 
was chosen because it requires fewer function evaluations than other methods for a given 
accuracy and does not require equal interval samples [Ref. 8]. CURVENEW, 
LOOPSCAT and HARLOOP read the coefficients from file OUTGLEG. The number 
of integration points per wavelength was varied to optimize convergence of LOOPSCAT 
and CURVENEW and is discussed in more detail in Chapter IV. The excitation vector 
[V] is calculated from equation (2.36) or (2.37) in subrouune CURVEW of 


SURVENEW. 


3. Execution Time 
Analysis of the nested DO loop structure of subroutine ZCURVED of 
program CURVENEW indicates that the total execution time of ZCURVED can be 


represented by, 


NN (3.5) 





Figure 6. Loop Geometry for Impedance Integrations 


Zh 


where N, is the number of curved segments (from equation (3.2)) N,, is the number of 
Gaussian integration constants per curved segment and a, is a constant. Execution time 
of the subroutines DECOMP and SOLVE, common to both CURVENEW and 


LOOPSCAT, can be represented similarly by [Ref. 9], 
1 = Vie 


where N is the number of segments. The excitation subroutine CURVEW execution time 


and field integrations are given by, 
= (NN Ga 


where again y and &, are constants. Assuming that the execution time of the rest of the 


program is negligible, the total execution time of CURVENEW is 


T, = T,,+T,+T;, = «,.N.N,.+¥No+0.N.N,. - (3.8) 


Cc 


A similar expression for LOOPSCAT uses the subscript 1, 
ee aN Noi +N; +,N/Noy =) 


Run times were recorded for various values of N., N, and N, using an IBM PC/AT with 
a math coprocessor and the coefficients for CURVENEW are found to be y=0.000156, 
a.=Q.0230, and ,=0.0222. The coefficients for LOOPSCAT are a,=0.0132 and 
£,=Q.0252. For the moment, assume that the number of Gaussian integration points per 
wavelength, N,, is held constant for both CURVENEW and LOOPSCAT. The number 


of integration points on a segment is, 
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Naa mv AL (3.10) 
and the number of segments 1s, 


2 
N. = ~ (Bain 





Similar expressions may be written for straight subsections. Combining equations (3.8), 
3.10, and 3.11 gives 
N-o 


c 3 
fe = 4n?re = +yN_+2t7r0N, ; (3.12) 


c 





iticeratio of I. to T, is given by, 


9 2 a72 s 
4n*roN,a.|N,+yN.+20r CN, 


Te (3.13) 


4n7rN,a,/N, +N; + 20790 ,N, 


Equation (3.13) will be used in Chapter IV to compare the execution times of 


CURVENEW and LOOPSCAT. 


IV. CALCULATED DATA FOR THE LOOP 


The convergence of the MM solutions for both the current and electric field for 
circular loops of various dimensions are presented for both linear and circular 
polarizations. The convergence is shown to depend on the segment size and number of 
integration points, as well as excitation conditions (incidence direction and polarization). 
Representative plots are presented within the chapter, and additional plots are given in 


Appendix B. 


A. CONVERGENCE OF HARLOOP 

The Fourier mode solution, HARLOOP, was tested for convergence with respect 
to incidence angle, number of modes, and number of integration constants to establish 
a baseline for comparison to CURVENEW and LOOPSCAT. HARLOOP was chosen as 
a baseline because the sinusoidal basis functions match the physical behavior of the 
current on the loop, and thus the current series converges rapidly. This is illustrated in 
Figure 7 for a 0.5 A radius loop with a plane wave incident at an angle of 40 degrees. 

Oscillation of the current as a function of ¢ becomes more rapid as @ is increased, 
because the phase of the incident field over the loop varies as sin © (equation (2.35)). 
For a @ polarized linear wave incident in the ¢=0 plane, the current is always zero at 
¢=0 and 180 degrees, where E' is cross-polarized with respect to the axis of the wire. 


For @ polarized incident waves, maxima occur at ¢=90 and 270 degrees, where E' is 
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x~107-2 Current Mognitude vs Loop angie 
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Figure 7. Convergence of the Current in the Complex Exponential Solution 


parallel to the axis of the wire. The overall current amplitude decreases for © 
approaching 90 degrees, as expected, since the loop’s projected area 1s small as viewed 
by the incident wave. For @¢ polarized incident waves, the minima occur at ¢=90 and 
270 degrees and maxima at ¢=0 and 180 degrees. The currents do not vanish for 0 
approaching 90 degrees because the loop is parallel to the ¢ polarized incident field. 
Circularly polarized incident waves give a constant magnitude current at normal 
incidence, and oscillations increase with 8. For 6=90 degrees, the circular and ¢ 
polarization responses are identical. 

HARLOOP is also found to be in agreement with measurements taken on the echo 
area of wire loops at normal incidence [Ref. 10]. The plot of Figure 8 gives the 
echo area (o/)’) versus rp for varying wire radius using HARLOOP. Measured data is 


indicated by the ’+’ sign. 


yes) 


B. CONVERGENCE OF THE CURRENT EXPANSION 

Having established HARLOOP as a baseline for comparison, convergence of the 
curved subsection program CURVENEW and the linear subsection program LOOPSCAT 
was evaluated. The plots of Figures 9 through 14 give the current on the loop as a 
function of loop angle, ¢, for varying ©, loop radius, and incident wave polarization. 
The number of integration points per wavelength, N,, is held constant at 320 in 
LOOPSCAT and CURVENEW. This number was chosen empirically to give a 
converged current within five to ten percent RMS. The RMS error is defined relative to 
the Fourier mode solution. The HARLOOP current is plotted with the solid line, those 
of CURVENEW are plotted with the "+" sign, and those of LOOPSCAT with the "o". 


The wire radius 1s 0.001 A for these calculations. 


Bockscotter Echo Areo va Loop Rodius 


Echo Area / wavelength~? 


wire rodiuax0.005 
inc. AnglesOceg 
— HARLOOPFP 
+ —- Meosured 
o,0OS | Ou Ore OrZo5 Gees fe ae 0.4 Ow a Ores 





Loop Rodius, wovelengths 


Figure 8. Backscatter Echo Area for a Loop with varying Radius at Normal Incidence 
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Figure 9. Magnitude of the Current on a 0.1 A Radius Loop, Normal 
Incidence, Circular Polarization (+ =curved; o =linear) 
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Figure 10. Phase of the Current on a 0.1 XA Radius Loop, Normal 
Incidence, Circular Polarization (+ =curved; o =linear) 
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Figure 11. Magnitude of the Current on a 0.1 X Radius Loop, Incidence 
Angle=40 deg, Circular Polarization (+ =curved; o =linear) 
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Figure 12. Phase of the Current on a 0.1] Radius Loop, Incidence 
Angle=40 deg., Circular Polarization (+ =curved; o =linear) 
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Figure 13. Magnitude of the Current on a 0.1 Radius Loop, Incidence 
Angle=40 deg., Theta Polarization (+ =curved; o =linear) 
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Figure 14. Phase of the Current on a 0.5 Radius Loop, Incidence 
Angle=40 deg., Theta Polarization (+ =curved; o =linear) 
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Representative plots of the normalized root mean squared error are given in Figures 15 
through 22. For set values of N and N,, CURVENEW converges faster than 
LOOPSCAT in most cases. The error difference 1s most pronounced for loop 
circumferences on the order of a wavelength or less. From the plots, for r,=0.1 A, 
CURVENEW converges to less than 10 percent error for segment sizes ranging from 
0.02 to 0.2 wavelengths (N,=32 to N,=3). LOOPSCAT converges to within 10 percent 
error for segment sizes less than approximately 0.06 A (N, > 10) but gives errors of 30 
to 40 percent for a segment size of 0.2 wavelengths. There is no improvement using 
curved subsections on larger loops for off-axis incidence waves, but the curved 
subsections give small improvements for large loops at normal incidence. This 1s 
expected in view of the behavior of the current on the loop. For linear polarization the 
current abruptly flips polarity from one side of the loop to the other. For circular 
polarization, the amplitude is constant, but the phase is linear. Both of these conditions 
can be represented accurately by a few triangles if the impedance and excitation integrals 
are evaluated precisely on the loop contour (see Figure 23). 

Plots of the magnitude of the backscattered E* versus incidence angle for varying 
segment size, N,, and ro are given in Figures 24 and 25. As with the currents, the 
backscattered field converges more rapidly for CURVENEW than LOOPSCAT in most 
cases, with the greatest difference for small radius loops and angles near the maximum 


values of | E*]. 
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Figure 15. Error in the Current Magnitude for a 0.1 > Radius Loop, 
Normal Incidence, Circular Polarization (+ =curved; o =linear) 
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Figure 16. Error in the Current Magnitude for a 0.1 A Radius Loop, 
Incidence Angle=40 deg., Circular Polarization (+ =curved; 0 =linear) 
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Figure 17. Error in the Current Magnitude for a 0.5 \ Radius Loop, 
Normal Incidence, Circular Polarization (+ =curved; o =linear) 
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Figure 18. Error in the Current Magnitude for a 0.5 X Radius Loop, 
Incidence Angle=40 deg., Circular Polarization (+ =curved; o =linear) 
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Figure 19. Error in the Current Magnitude for a 0.1 Radius Loop, 
Normal Incidence, Theta Polarization (+ =curved; 0 =linear) 
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Figure 20. Error in the Current Magnitude for a 0.1 \ Radius Loop, 
Angle of Incidence =40 deg., Theta Polarization (+ =curved; o =linear) 
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Figure 21. Error in the Current Magnitude for a 0.5 \ Radius Loop, 
Normal Incidence, Theta Polarization (+ =curved; o =linear) 
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Figure 22. Error in the Current Magnitude for a 0.5 X Radius Loop, 
Incidence Angle=40 deg., Theta Polarization (+ =curved; o =linear) 
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Figure 23. A Representation of a Sinusoidal Current with Two Basis Functions (top) and 
a Constant Current as a Superposition of Triangles (bottom) 
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C. COMPARISON OF EXECUTION TIME AND MEMORY REQUIREMENTS 

The average run time for a given loop radius and convergence error is greater for 
CURVENEW than for LOOPSCAT due to the N,’ dependence in equation (3.13) and 
relative magnitudes of the coefficients a and y. The plot of equation (3.13) in Figure 26 


illustrates the ratio of run times of CURVENEW and HARLOOP versus N, for N.=4, 


Ratio of Run Times vs. N! 
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Figure 26. Ratio of Execution Times for CURVENEW and LOOPSCAT for Fixed 
Curved Segment Lengths and with Varying Linear Segment Lengths 


8, 32 and N,=320 for a 0.1 A loop. The run time of CURVENEW is less than 
HARLOOP for T,/T, < 1. From the plot, the "break even" points are approximately 
N,=115, 90, 52 for N.=4, 8, 32. For N, less than these values, the integration 
subroutine ZCURVED is the determining factor in the run time; for N, greater than these 
values, Gausssian elimination subroutines DECOMP and SOLVE are the determining 


factors. The increased number of integrations per segment completely offset the time 


oF 


savings of a reduced [Z] matrix in CURVENEW. As mentioned in Chapter III, a delta 
function approximation for the outer integration of the impedance integral of equation 
(2.28) was investigated. This reduces the exponent of N, to one in equation (3.13), but 
many more segments are required for a given accuracy. A more efficient integration 
scheme using a large number of integration points in the vicinity of 6=¢' and fewer 
integrations elsewhere may reduce the execution time. 

The savings in computer memory is significant for CURVENEW, since the number 
of matrix elements is on the order of N’. For a given accuracy for a 0.1 loop, the ratio 
of the number of elements required for curved and linear subsections, (N./N,)’, is on the 
order of 0.1 to 0.2. This is a reduction of 80 to 90 percent. Although the memory 
requirements for the small loops considered here are not prohibitive for piecewise linear 
segments, greater memory savings will be realized for larger geometries comprised of 


many small features. 
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V. CONCLUSIONS 


The use of conformal subdomain basis functions (curved subsections) to represent 
the current on a thin curved wire was investigated by solving the thin wire electric field 
integral equation using the method of moments. A solution using triangular basis 
functions was computer coded in FORTRAN and validated by comparing it to measured 
data and the results of two other method of moments solutions (LOOPSCAT and 
HARLOOP). The effect of varying loop radius, segment size, number of integration 
points and incident wave parameters on the accuracy and rate of convergence of the 
Current expansion and backscattered field was investigated. 

For small loops with circumferences on the order of a wavelength, the number of 
segments required to converge to a given accuracy with the curved segments was as 
small as 20 percent of the number of linear segments required to converge to the same 
accuracy (see Table 2). From computed data, it was determined that the greatest 
reduction in the number of unknowns for curved subsections occurs for geometries where 
the current amplitude variations over the surface are small and the phase variations are 
small or linear. As mentioned in Chapter I, the general rule of thumb for the length of 
one segment is 0.05 A to 0.1 A, which corresponds to a phase variation of 20 to 40 
degrees. This restriction in phase variation is the driving factor when choosing the 
segment size for geometrical features having radii of curvature of approximately /2 or 


larger. A piecewise linear segment is small enough to represent the geometry accurately 
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in this case. The phase restriction applies to curved segments as well, but the curved 
segments conform exactly to the wire, and hence for small loops there is no sacrifice in 
geometrical accuracy by choosing segment sizes of approximately 0.05 A or 20 electrical 
degrees. As the loop becomes larger, or the wavelength becomes smaller, the curved and 
straight subsection solutions become equivalent. The greatest advantage in using curved 
subsections to reduce the number of segments 1s for electrically small structures where 
small linear segments are required simply to reproduce the wire shape. 

Although the number of segments was greatly reduced using conformal subsections, 
the execution time was increased due to the increased number of integration points per 
segment required for acceptable accuracy. To reduce the integration time, it is suggested 
that the number of integration points per wavelength be varied from a large number when 
evaluating the self term, to fewer points away from the self term. For certain geometries, 
symmetry could also be used to reduce the integration time. 

To avoid singularities, the MM testing procedure was performed along the axis of 
the wire and the current constrained to the surface of the wire. A delta function 
approximation in the impedance integrations was found to reduce the time required to 
compute the impedance matrix, but as expected, required more segments for a given 
convergence accuracy. 

A disadvantage in the formulation of CURVENEW is that an analytic expression 
for the curve is needed to perform the integrations for the impedance matrix and the 
excitation vector. The expressions will change each time a new curve is analyzed, and 


consequently, considerable effort will be required to modify the code every time a change 
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is made in the geometry. Programs that use linear segments are more flexible because 
they require only the coordinates of the points along the wire axis to generate the 
integration points. 

The next logical step in testing the effectiveness of conformal subdomain basis 
functions is to formulate the solution for an equiangular spiral wire. The equiangular 
spiral is used in broadband antennas and has a simple mathematical form. For 
geometrical accuracy, the segment size in the piecewise linear formulation will be much 


smaller than 20 electrical degrees near the center of the spiral. Equal length conformal 


TABLE 2. COMPARISON OF CURVENEW AND LOOPSCAT 


a a a 







Curved Linear | 
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Number of Segments for 10% 
RMS Error, Off Axis Incidence 


Execution Time’, 10% RMS 
Error, Normal Incidence 


Execution Time’, 10% RMS 
Error, Off Axis Incidence 


[Z] Matrix Size, 10 % RMS 
Error, Normal Incidence 


| [Z} Matnx Size, 10 % RMS 
Error, Off Axis Incidence 






















* Execution time measured with an IBM PC/AT. 


4] 


segments may be used along the spiral arms, and it is anticipated that the number of 


required segments will be substantially reduced. 
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APPENDIX A 


COMPUTER CODES 


C MAIN PROGRAM *CURVNEW.FOR* 

C **** MORE NUMERIGAELY EFFICIENT SEAN CURT SUE = —- 

C PLANE WAVE SCATTERING FROM A CIRCULAR LOOP IN THE Z PLANE. 
C METHOD OF MOMENTS WITH CURVED SUBSECTIONS 

C **** RADAR CROSS SECTION CALCULATION ==" 


@ 
COMPLEX Z(25000),B(500),C(500),BP(500),CP(500),EX,EC 
COMPLEX ET,EP,UC,RW(500),U0 
DIMENSION ECP(400),IPS(500), ANG(400), EDP(400) 
DIMENSION XG(128),AG(128),PH(500), DEL(500),PA(500) 
DIMENSION ECV(500), EX V(500), PHC(500),PHX(500) 
DATA PI/3.14159265/ 

DATA IPRINT/1/,ITEST/1/ 
RAD=PI/180. 

ECX=0. 

BK =2.*P] 

ETA =377. 

U0=(0.,0.) 
UC=(0.,-1.)*ETA*BK/4./PI 
CONST = 16.*PI**3 
PHIRD=0. 

C 

C READ INPUT AND PROGRAM CONTROL PARAMETERS 

C 
OPEN(1,FILE='PARAMLST. DAT’) 

READ(1,*) ANGLE, DT,START, STOP, AW,RO,SEG, HARR,GCONST,XLOC, XIPOL 
LOC=INT(XLOC) 

IPOL =INT(XIPOL) 

CLOSE(1) 

eC 

C READ GAUSSIAN CONSTANTS 

c 
OPEN(2,FILE=’OUTGLEG’) 

READ(2,*) NG 
DO 1 I=1,NG 
READ(2,*) XG(I), AG(I) 
1 CONTINUE 
C 


C OUTCURV IS THE FILE THAT THE SCATTERED FIELD DATA IS WRITTEN TO 
C ICURVOUT IS THE FILE THAT THE LOOP CURRENT DATA IS WRITTEN TO 
€ 
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43 
44 
45 
46 
47 
48 
49 
50 
51 
52 
53 
54 
55 
56 
57 
58 
59 
60 
61 
62 
63 
64 
65 
66 
67 
68 
69 
70 
71 
72 
73 
74 
75 
16 
77 
78 
79 
80 
81 
82 
83 
84 
85 
86 
87 
88 
89 
90 
91 
92 


OPEN(8,.FILE=’OURCURV.DAT’) 
OPEN(7,FILE=*ICURVOUT.DAT’) 
@ 
C GENERATE THE LOOP POINTS 
€ 
C 
C CHOOSE THE NUMBER OF POINTS BASED ON THE VALUE OF SEG 
eG 
C AW=AW*RO 
DPHI=SEG/RO 
NP=INT(2.*PI/DPHI) +1 
DPHI=360./FLOAT(NP) 
PH(1)=0. 
DO 10 [=2,NP+1 
PH(I)=FLOAT(I-1)*DPHI*RAD 
DEL(I-1)=(PH(I)-PH(I-1)) 
PA(I-1)=(PH(1) + PH(I-1))/2. 
10 CONTINUE 
NP=NP+2 
C 
C OVERLAP THE ENDS SO THAT CURRENT WILL BE CONTINUOUS ON THE LOOP 
E 
PH(NP)=BK + PH(2) 
DEL(NP-1)=DEL(1) 
PA(NP-1)=BK+PA(1) 
MT=NP-2 
DO 52 I=1,NP 
XHB=R0*COS(PH(1)) 
YHB=RO*SIN(PH(I)) 
52. CONTINUE 
WRITE(6,*) "GEOMETRY DEFINED’ 
IF(ITEST.EQ.0) GO TO 98 
C 
C DEFINE DIMENSIONS OF THE IMPEDANCE MATRIX BLOCKS 
6 
WRITE(6,*) ‘NP,MT=",NP,MT 
G 
C COMPUTE IMPEDANCE MATRIX ELEMENTS 
C 
CALL ZCURVED(NP,RO,PH,DEL,PA,NG,XG,AG,AW,Z,ZMAX) 
DO 11 I=1,MT 
CZ=CABS(Z(I)) 
AZ=ATAN2(AIMAG(Z(I)), REAL(Z(1)) + 1.E-8)/RAD 
11 CONTINUE 
WRITE(6,*) *WIRE IMPEDANCE COMPUTED’ 
C 
C PERFORM LU DECOMPOSITION 
C 
CALL DECOMP(MT, IPS,Z) 
WRITE(6,*) ’Z DECOMPOSED’ 
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93 
94 
95 
96 
97 
98 
99 
100 
10] 
102 
103 
104 
105 
106 
107 
108 
109 
110 
11] 
112 
ts 
114 
Hie 
116 
117 
118 
eee 
120 
ed 
122 
123 
124 
125 
126 
i eay, 
128 
9 
130 
ea 
132 
3 
134 
135 
136 
be 
138 
139 
140 
14] 
142 


Cc 
C BEGIN FIELD CALCULATIONS 
C 
98 PHRO=PHIRD*RAD 
IT=INT((STOP-START)/DT) + 1 
WRITE(7,*) IT,MT,0,0 
DO 500 I=1,IT 
THETA =FLOAT(I-1)*DT +START 
THR=THETA*RAD 
PHR=PHRO 
IF(THETA.LT.180.) GO TO 99 
THR =(360.-THETA)*RAD 
PHR=PHRO+PI 
99 CONTINUE 


ET=UO0 
EP=U0 
c 
C COMPUTE THE EXCITATION VECTOR 
C 
CALL CURVEW(NP,RO,PH,DEL,PA,NG,XG,AG,THR,PHR,RW) 
lf (LOC EQO0) THEN 
c 


C CIRCULAR POLARIZATION IF LOC=1 ELSE LINEAR 
: IF(IPOL.EQ.1) THEN 

: THETA POLARIZED INCIDENT WAVE (JPOL=1) 

' DO 101 L=1,MT 


101 B(L)=RW(L) 
ELSE 
G 
C PHI POLARIZED INCIDENT WAVE (IPOL =2) 
c 
ENDIF 


IF(ITEST.EQ.0) GO TO 9998 
Cc 
C PERFORM GAUSSIAN ELIMINATION TO DETERMINE [C} 
C 
CALL SOLVE(MT, IPS,Z,B,C) 
DO 210 L=1,MT 
WRITE(7,*) L,L*DPHI,CABS(C(L)), ATAN2(REAL(C(L)), 
* AIMAG(C(L)))/RAD 
ET=ET+RW(L)*C(L) 
210 EP=EP+RW(L+MT)*C(L) 


EESE 
Cc 
C THETA POLARIZED INCIDENT WAVE 
C 


DOZ7 iE ihm 
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143 
144 
145 
146 
147 
148 
149 
150 
151 
152 
153 
154 
155 
156 
157 
158 
159 
160 
161 
162 
163 
164 
165 
166 
167 
168 
169 
170 
17] 
ee 
173 
174 
175 
176 
177 
178 
leo 
180 
18] 
182 
183 
184 
185 
186 
187 
188 
189 
190 
Eo 
os 


Zo 


S 


B(L)=RW(L) 


C PHIEPOLARIZED INCIDENT WAVE 
C PHASE SHIFT FOR CPUS Pi/2. 


C 


222 


Zeu 


DO 222 L=1,MT 
BP(L)=RW(L+ MT)*CEXP(CMPLX(0.,PI/2.)) 

CALL SOLVE(MT, IPS,Z,B,C) 

CALL SOLVE(MT, IPS,Z,BP,CP) 

DO 220 L=1,MT 
WRITE(7,*) L,L*DPHI,CABS(C(L) + CP(L)), ATAN2(REAL(C(L) + CP(L)), 
AIMAG(C(L) + CP(L)))/RAD 
ET=ET+RW(L)*C(L) + RW(L)*CP(L) 

EP=EP+RW(L+MT)*C(L) +RW(L+MT)*CP(L) 

ENDIF 

EC=ET*UC 

EX =EP*UC 

ANG(I)=THETA 

ECV(1)=CABS(EC) 

EX V(I) = CABS(EX) 

ECR=REAL(EC) 

ECI= AIMAG(EC) 

EXR=REAL(EX) 

EXI = AIMAG(EX) 

PHC(1) = ATAN2(ECI,ECR + 1.E-20)/RAD 

PHX(I)=ATAN2(EXI1,EXR + 1.E-20)/RAD 

ECX = AMAX I(ECX, ECV(1), EX V(J)) 


500 CONTINUE 


WRITE(6,*) “EMAX=’,ECX 
DO 600 I=1,IT 


ECV(I)=AMAXI(ECV(I), 1.E-10) 
EXV(I)=AMAXI(EXV(I), 1.E-10) 
ECP(1) = (ECV(I)/ECX)**2 
EDP(I)=(EXV(1)/ECX)**2 
ECP(1)=AMAX1(ECP(1),.00001) 
EDP(I)= AMAX1(EDP(1),.00001) 
ECP(I) = 10.*ALOG10(ECP(1)) 
EDP(I)=10.*ALOG 10(EDP(1)) 


600 CONTINUE 


SIGMA =(ECX**2)*CABS(UC)/(2.*BK) 
SIGDB= 10.*ALOG10(SIGMA) 
WRITE(6,*) "BACKSCATTER CROSS-SECTION, IN DB=’,SIGMA,SIGDB 


208 FORMAT(/,5X,’SIGMA/WAVL SQ=’,E]5S.4, 


Bao. G IN DB=’,F8.4) 
DO 9000 L=1,IT 


WRITE(8,*) ANG(L),ECV(L) 


9000 CONTINUE 
9998 35 TOP 


END 
SUBROUTINE SOLVE(N,IPS,UL,B,X) 
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220 


bo 


e 
Cc 


C SUBROUTINE TO SOLVE SYSTEM OF EQUATIONS WITH COMPLEX 


COEFFICIENTS: 
C CALL ’DECOMP’ FIRST. (FROM MAUTZ AND HARRINGTON) 


S 


Cc 


COMPLEX UL(50000),B(200),X(200),SUM 
DIMENSION IPS(500) 
NP1=N+1 

IP=IPS(1) 

X(1)=B(IP) 

DO 2 I=2,N 

IP=IPS(I) 

IPB=IP 

IM1=I-1 

SUN =(0.,0,) 

DO 1 J= Ini 
SUM=SUM + UL(IP)*X(J) 


PIP—IesoN 
2 X(1)=B(IPB)-SUM 


K2=N*(N-1) 
IP=IPS(N)+K2 
X(N)=X(N)/UL(IP) 
DO 4 IBACK=2,N 
I=NPI-IBACK 
K2=K2-N 
IPI=IPS()+K2 
Wei 
SUM=(0.,0.) 
Jeeta 

DO 37 —IP iN 
IP=IP+N 


3 SUM=SUM + UL(IP)*X(J) 
4 X(1I)=(X(1)-SUM)/UL(PI) 


RETURN 
END 
SUBROUTINE DECOMP(N,IPS,UL) 


C SUBROUTINE TO DECOMPOSE SYSTEM OF EQUATIONS. 
C FROM MAUTZ AND HARRINGTON. 


Cc 


COMPLEX UL(50000),PIVOT,EM 

DIMENSION SCL(200),IPS(200) 

DO 5 I=1,N 

IP Stij—1 

RN=0. 

Visi 

DO 2 J=1,N 

ULM=ABS(REAL(UL(J1))) + ABS(AIMAG(UL(J1))) 
Ji=jieeN 
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243 PP CRING OEM ) 1.2.2 


244 1 RN=ULM 

245 2 CONTINUE 

246 SCL(I) =1./RN 

247 5 CONTINUE 

248 NMI=N-1 

249 K2=0 

250 DO 17 K=1,NM1 

251 BIG=0. 

252 DO 11 I=K,N 

253 IP=IPS(I) 

254 IPK=IP+K2 

255 SIZE =(ABS(REAL(UL(IPK))) + ABS(AIMAG(UL(IPK))))*SCL(IP) 
256 IF(SIZE-BIG) 11,11,10 

257 10 BIG=SIZE 

258 IPV=I 

259 11 CONTINUE 

260 IF(IPV-K) 14,15,14 

261 14 J=IPS(K) 

262 IPS(K) =IPS(IPV) 

263 IPS(IPV)=J 

264 15 KPP=IPS(K)+K2 

265 PIVOT=UL(KPP) 

266 KPI=K+1 

267 DO 16 I=KP1,N 

268 KP=KPP 

269 IP=IPS(I)+K2 

270 EM=-UL(IP)/PIVOT 

271 18 UL(IP)=-EM 

72 DO 16 J=KP1,N 

oe IP=IP+N 

274 KP=KP+N 

275 UL(IP)= UL(IP) + EM*UL(KP) 

276 16 CONTINUE 

oT K2=K2+N 

278 17 CONTINUE 

279 RETURN 

280 END 

281 SUBROUTINE ZCURVED(NP,RO,PH,DEL,PA,NG,XG,AG,A,Z,ZMAX) 
282 C 

283 C IMPEDANCE ELEMENTS FOR CURVED BASIS FUNCTIONS. 
284 C SPECIFICALLY DERIVED FOR A CIRCULAR LOOP -- NOT A GENERAL CASE. 
285 G 

286 COMPLEX CEXP,Z(50000),CON,CMPLX,SUMA,SUMB 
287 COMPLEX U0,SUM1,SUM2,SUM3,SUM4,EXP,ZT(500) 
288 DIMENSION DEL(500),PH(500),XG(128), AG(128),PA(500) 
289 C OPEN(2,FILE=’ZCURV.DAT’) 

290 ETA =377. 

291 ZMAX=0. 

292 PI=3.14159 
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203 
294 
295 
296 
e77 
298 
292 
300 
301 
302 
303 
304 
305 
306 
307 
308 
309 
310 
S11 
312 
313 
314 
315 
316 
317 
318 
219 
320 
321 
322 
a03 
324 
325 
326 
a2 / 
328 
a2 
330 
33] 
332 
333 
334 
335 
336 
a3/ 
338 
339 
340 
341 
342 


BK =2.*PI 
BK2=BK**2 
U0=(0.,0.) 
CON =(0.,1.)*BK*ETA/(4.*PI)*RO**2 
NT =NP-2 
C 
C COMPUTE Z(1,LQ) = ZT(LQ) 
@ 
KQ=1 
P1=DEL(KQ)/2. 
P2=PA(KQ) 
P3=DEL(KQ +1)/2. 
P4=PA(KQ+1) 
C 
C DO THE L Loop 
G 
DO 600 LQ=1,NT 
PP1=DEL(LQ)/2. 
PP2=PA(LQ) 
PP3 =DEL(LQ+1)/2. 
PP4=PA(LQ+1) 
@ 
C DO THE PHI INTEGRATION 
C *** FIRST PART FROM PHI(K) TO PHI(K +1) 
C PHI PRIMED INTEGRATION FOR THE POSITIVE SLOPE OF LQ 
C 
SUMA=U0 
PHA =PA(KQ) 
DO 100 I=1,NG 
PHI=P1*XG(I) + P2 
TK =(PHI-PH(KQ))/DEL(KQ) 
TKP=1./DEL(KQ)/RO 
SUM1=U0 
DO 90 J=1,NG 
PHIP=PP1*XG(J) + PP2 
TL =(PHIP-PH(LQ))/DEL(LQ) 
TLP=1./DEL(LQ)/RO 
CC=COS(PHI-PHIP) 
€ 
C COMPUTE THE MAGNITUDE OF R. NOTE THAT R IS COMPUTED FROM THE 
C WIRE AXIS TO THE SURFACE OF THE WIRE 
eC 
RR=RO*SQRT(4.*(SIN((PHI-PHIP)/2.))**2 + (A/RO)**2) 
EXP=CEXP(CMPLX(0.,-BK*RR))/RR 
SUM1=SUM1+AG(J)*EXP*(TK*TL*CC-TKP*TLP/BK2) 
90 CONTINUE 
SUM1=SUM1*PP1 
C 
C PHI PRIMED INTEGRATION FOR THE NEGATIVE SLOPE OF LQ 
c 
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343 
344 
345 
346 
347 
393 
394 
395 
396 
397 
398 
399 
400 
401 
402 
403 
404 
405 
406 
407 
408 
409 
410 
41] 
412 
413 
414 
415 
416 
417 
418 
419 
420 
42] 
422 
423 
424 
425 
426 
427 
428 
429 
430 
431 
432 
433 
434 
435 
436 
437 


SUM2=U0 
DO 80 J=1,NG 
PHIP =PP3*XG(J) + PP4 
TL =1.-(PHIP-PH(LQ + 1))/DEL(LQ +1) 
TLP=-1./DEL(LQ+1)/RO 
CC=COS(PHI-PHIP) 
E 
C COMPUTE THE MAGNITUDE OF R. NOTE THAT R IS COMPUTED FROM THE 
C WIRE AXIS TO THE SURFACE OF THE WIRE 
C 
RR=RO*SQRT(4. *(SIN((PHI-PHIP)/2.))**2 + (A/RO)**2) 
EXP=CEXP(CMPLX(0.,-BK*RR))/RR 
SUM2=SUM2+ AG(J)*EXP*(TK*TL*CC-TKP*TLP/BK2) 
80 CONTINUE 
SUM2=SUM2*PP3 
SUMA=SUMA +(SUM1+SUM2)*AG(I) 
100 CONTINUE 
SUMA=SUMA#*P1 
C 
C *** SECOND PART FROM PHI(K+1) TO PHI(K +2) 
C PHI PRIMED INTEGRATION FOR THE POSITIVE SLOPE OF LQ 
C 
SUMB=U0 
PHA=PA(KQ +1) 
DO 101 I=1,NG 
PHI=P3*XG(I)+P4 
TK =1.-(PHI-PH(KQ+ 1))/DEL(KQ+1) 
TKP=-1./DEL(KQ+1)/RO 
SUM3=U0 
DO 91 J=1.NG 
PHIP=PP1*XG(J)+PP2 
TL=(PHIP-PH(LQ))/DEL(LQ) 
TLP=1./DEL(LQ)/RO 
CC=COS(PHI-PHIP) 
‘es 
C COMPUTE THE MAGNITUDE OF R. NOTE THAT R IS COMPUTED FROM THE 
C WIRE AXIS TO THE SURFACE OF THE WIRE 
G 
RR =RO*SQRT(4. *(SIN((PHI-PHIP)/2.))**2 + (A/RO)**2) 
EXP=CEXP(CMPLX(0.,-BK*RR))/RR 
SUM3 =SUM3 + AG(J)*EXP*(TK*TL*CC-TKP*TLP/BK2) 
9] CONTINUE 
SUM3=SUM3*PP1 
C 
C PHI PRIMED INTEGRATION FOR THE NEGATIVE SLOPE OF LQ 
E 
SUM4=U0 
DO 81 J=1,NG 
PHIP=PP3*XG(J)+PP4 
TL=1.-(PHIP-PH(LQ + 1))/DEL(LQ +1) 


51 


493 
494 
495 
496 
497 
498 
499 
500 
501 
502 
503 
504 
505 
506 
507 
508 
509 
10 
Sid 
a2 
513 
514 
515 
516 
2) 
518 
519 
520 
a2] 
Dee 
23 
524 
25 
526 
a7 
528 
529 
530 
531 
32 
533 
534 
235 
536 
537 
538 
339 
540 
541 
542 


TLP=-1./DEL(LQ+1)/RO 
CC=COS(PHI-PHIP) 
c 
C COMPUTE THE MAGNITUDE OF R. NOTE THAT R IS COMPUTED FROM THE 
C WIRE AXIS TO THE SURFACE OF THE WIRE 
G 
RR=RO*SQRT(4.*(SIN((PHI-PHIP)/2.))**2 + (A/RO)**2) 
EXP =CEXP(CMPLX(0.,-BK*RR))/RR 
SUM4=SUM4+AG(J)*EXP*(TK*TL*CC-TKP*TLP/BK2) 
8] CONTINUE 
SUM4=SUM4*PP3 
SUMB=SUMB +(SUM3 +SUM4)*AG(I) 
101 CONTINUE 
SUMB=SUMB*P3 
ZT(LQ)=CON*(SUMA + SUMB) 
ZMAX = AMAX1(ZMAX, CABS(ZT(LQ))) 
600 CONTINUE 
ZT(NT) =ZT(2) 
G 
C FILL THE ENTIRE Z MATRIX USING SYMMETRY PROPERTIES 
C [Z] IS A SYMMETRICAL TOEPLITZ MATRIX 
C ROW INDEX, I; COL INDEX, J 
Cc 
DO 10 I=1,NT 
DO 10 J=1,NT 
K =(I-1)*NT+J 
Z(K)=U0 
UJ =IABS(I-J) 
IF(IJ.GT.NT) GO TO 10 
Ui=U4+1 
Z(K)=ZT(IJ1) 
10 CONTINUE 
RETURN 
END 
SUBROUTINE CURVEW(NP,RO,PH,DEL,PA,NG,XG,AG,THR,PHR,R) 
e 
C PLANE WAVE EXCITATION VECTOR ELEMENTS FOR A LOOP USING CURVED 
C BASIS FUNCTIONS. INCIDENCE DIRECTION IS (THR,PHR). THE WIRE LIES 
C IN THE X-Y PLANE (Z=0). 
€ 
COMPLEX UO0,R(500),SUM1,SUM2,SUM3,SUM4,CEXP, FF 
COMPLEX GG,CMPLX,SUMT,SUMP 
DIMENSION PH(500),DEL(500),PA(500),XG(128), AG(128) 
MT =NP-2 
U0=(0.,0.) 
PI =3.14159 
BK =2.*PI 
ST=SIN(THR) 
CT=COS(THR) 
DO 50 IP=1,MT 
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543 
544 
545 
546 
547 
548 
549 
550 
a2) 
Joe 
523 
554 
325) 
556 
a3) 
558 
Sade 
560 
561 
562 
563 
564 
565 
566 
567 
568 
507 
570 
1 
oy 4 
Se 
574 
575 
576 


SUM1=U0 
SUM2=U0 
SUM3=U0 
SUM4=U0 
P1=DEL(IP)/2. 
P2=PA(IP) 
P3=DEL(IP+1)/2. 
P4=PA(IP +1) 
DO 20 I=1,NG 
PHI=P1*XG(I)+P2 
CC=COS(PHR-PHI) 
SS=SIN(PHR-PHI)*CT 
FF = AG(I)*(PHI-PH(IP))/DEL(IP)*CEXP(CMPLX(0.,BK*RO*ST*CC)) 
SUM1=SUM1+CC*FF 
SUM2=SUM2+SS*FF 
PHI=P3*XG(I)+P4 
CC=COS(PHR-PHI) 
SS =SIN(PHR-PHI)*CT 
GG = AG(I)*(1.-(PHI-PH(IP + 1))/DEL(IP + 1))*CEXP(CMPLX(0., 
* BK*RO*ST*CC)) 
SUM3=SUM3+CC*GG 
SUM4=SUM4+SS*GG 
20 CONTINUE 
SUMP=SUM1*P1 + SUM3*P3 
SUMT =SUM2*P1 +SUM4*P3 
‘6 
C R-WIRE-THETA IN R(IP) AND R-WIRE-PHI IN R(IP+MT) 
re 
R(IP)=SUMT*RO 
R(IP+MT)=SUMP*RO 
50 CONTINUE 
RETURN 
END 
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C MAIN PROGRAM *LOOP.FOR* 
C PLANE WAVE SCATTERING FROM A CIRCULAR LOOP IN THE Z PLANE. 
C **** RADAR CROSS SECTION CALCULATION ****# 
Cc 
COMPLEX Z(15000),B(500),C(500),BP(500),CP(500),U 
COMPLEX ET,EP,UC,RW(500),U0 
DIMENSION ECP(400),IPS(500), ANG(400), EDP(400) 
DIMENSION ZH(200),XT(128),AT(128), XH(200), YH(200) 
DIMENSION ECV(400),EXV(400),PHC(400),PHX(400) 
DATA PI/3.14159265/ 
DATA IPRINT/1/ 
C 
C READ INPUT AND PROGRAM CONTROL PARAMETERS 
C 
OPEN(1,FILE=’PARAMLST.DAT’) 
READ(1,*)ANGLE,DT,START,STOP, AW,RB,SEG,HARR,GCONST,XLOC, XIPOL 
LOC =INT(XLOC) 
IPOL=INT(XIPOL) 
CLOSE(1) 
G 
C READ GAUSSIAN CONSTANTS 
C 
OPEN(2,FILE=’OUTGLEG’) 
READ(2,*) NT 
DO 2 J=1,NT 
READ(2,*) XT(I),AT(I) 
2 CONTINUE 
RAD =PI/180. 
ECX =0. 
BK =2.*P] 
ETA=377. 
U=(0.,1.) 
U0=(0.,0.) 
UC=-U*ETA*BK/4./PI 
CONST = 16.*PI**3 
NT2=NT/2 
E 
C OUTLOOP IS THE FILE THAT THE SCATTERED FIELD DATA IS WRITTEN TO 
C ISTOUT IS THE FILE THAT THE LOOP CURRENT DATA IS WRITTEN TO 
C 
OPEN(8,FILE='OUTLOOP.DAT’) 
OPEN(7,FILE=*ISTOUT. DAT’) 
DPHI=SEG/RB 
NP=INT(2.*PI/DPHI) +1 
DPHI =360./FLOAT(NP) 
WRITE(6,*) ‘DPHI,NP=',DPHI,NP 
e 
C GENERATE THE LOOP POINTS. MULTIPLY ALL QUANTITIES BY BK (=2*PI) 
C 
C 
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52 
53 
54 
55 
56 
57 
58 
59 
60 
61 
62 
63 
64 
65 
66 
67 
68 
69 
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1 
74 
75 
16 
77 
78 
79 
80 
8] 
82 
83 
84 
85 
86 
87 
88 
89 
90 
91 
92 
93 
94 
95 
96 
97 
98 
99 
100 


C CHOOSE THE NUMBER OF POINTS BASED ON THE VALUE OF SEG 
€ 
AK =AW*BK 
DO 10 I=1,NP+1 
PP =FLOAT(I-1)*DPHI*RAD 
XH(1) =RB*COS(PP)*BK 
YH(1)=RB*SIN(PP)*BK 
ZH(1)=0. 
10 CONTINUE 
NP=NP+2 
c 
C OVERLAP THE ENDS SO THAT THE CURRENT IS CONTINUOUS 
GC 
XH(NP)=XH(2) 
YH(NP)= YH(2) 
ZH(NP) = ZH(2) 
DO 52 I1=1,NP 
YHB=YH(l)/BK 
XHB=XH(l)/BK 
DEL=0. 
IF(I.NE.1) THEN 
DXX = XHB-XH(I-1)/BK 
DYY = YHB-YH(I-1)/BK 
DEL=SQRT(DXX**2+DYY**2) 
ENDIF 
52. CONTINUE 
¢ 
C DEFINE DIMENSIONS OF THE IMPEDANCE MATRIX BLOCKS 
C 


MT =NP-2 

WRITE(6,*) "MT=',MT 
8: 
CCOMP UTE IMPEDANCE MATRIXCELEMENTS 
e 


CALL ZMATWW(1.1,NP,RB,XH,YH,ZH,NT.XT,AT,AK,Z) 
WRITE(6,*) *Z COMPUTED’ 
IF(IPRINT.EQ.0) THEN 
DO 1010I=1,MT 
CZ = CABS(Z(1)) 
AZ=ATAN2(AIMAG(Z(I)), REAL(Z(I)) + 1.E-8)/RAD 
WRITE(6,*) "I, Z=’,1,Z(1),CZ,CZ/ZMAX,AZ 
1010 CONTINUE 


ENDIF 
G 
C PERFORM LU DECOMPOSITION 
C 
CALL DECOMP(MT, IPS,Z) 
WRITE(6,*) °Z DECOMPOSED’ 
€ 


€ BEGIN TIERED CALCULATIONS. PHI FOR PATTERN CUT (DEGREES)=PHID 
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101 
102 
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104 
105 
106 
107 
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109 
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143 
144 
145 
146 
147 
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149 
150 


PHID=0. 
PHRO=PHID*RAD 
IT=INT((STOP-START)/DT) +1 
WRITE(7,*) IT,MT,0,0 
DO 500 I=1,IT 
THETA =FLOAT(I-1)*DT+START 
THR=THETA*RAD 
PHR=PHRO 
IF(THETA.LT.180.) GO TO 99 
THR = (360.-THETA)*RAD 
PHR=PHRO+PI 
99 CONTINUE 


ET=U0 
EP=U0 
C 
C COMPUTE THE EXCITATION VECTOR 
C 
CALL PLANEW(NP,XH, YH,ZH, THR,PHR,RW) 
IF (EOG7 EO 0) THEN 
6 


C CIRCULAR POLARIZATION IF LOC=1 ELSE LINEAR 
€ 
IF(IPOL.EQ.1) THEN 
€ 
C THETA POLARIZED INCIDENT WAVE (IPOL=1) 
€ 
DO 101 L=1,MT 
101 B(L)=RW(L) 
ELSE 
C 
C PHI POLARIZED INCIDENT WAVE (IPOL=2) 
C 
DO 102 L=1,MT 
102 B(L)=RW(L+ MT) 
ENDIF 
C 
C PERFORM GAUSSIAN ELIMINATION TO DETERMINE [C}] 
G 
CALL SOLVE(MT,IPS,Z,B,C) 
DO 210 L=1,MT 
WRITE(7,*) L,L*DPHI,CABS(C(L)), ATAN2(REAL(C(L)), 
* AIMAG(C(L)))/RAD 
ET =ET+(RW(L)/BK)*C(L) 
210 EP=EP+(RW(L+MT)/BK)*C(L) 


BESe 
GC 
C THETA PLOARIZED INCIDENT WAVE 
C 


DO 221 L=1,MT 
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15] 
2 
bd3 
154 
155 
156 
157 
158 
159 
160 
161 
162 
163 
164 
165 
166 
167 
168 
169 
170 
171 
2 
173 
174 
175 
176 
heey 
178 
179 
180 
181 
182 
183 
184 
185 
186 
187 
188 
189 
190 
19] 
192 
193 
194 
195 
196 
197 
198 
199 
200 


221 + B(L)=RW(L) 
c 
C PHI POLARIZED INCIDENT WAVE 
C PHASE SHIFT FOR CP IS PI/2. 
CG 
DO 222 L=1,MT 
222. ~=BP(L)=RW(L+MT)*CEXP(CMPLX(0. ,PI/2.)) 
CALL SOLVE(MT,IPS,Z,B,C) 
CALL SOLVE(MT, IPS,Z,BP,CP) 
DO 220 L=1,MT 
WRITE(7,*) L,L*DPHI, CABS(C(L) + CP(L)), ATAN2(REAL(C(L) + CP(L)), 
*  AIMAG(C(L) +CP(L)))/RAD 
ET=ET+RW(L)*C(L)+RW(L)*CP(L) 
220 EP=EP+RW(L+MT)*C(L)+RW(L+MT)*CP(L) 


ENDIF 

ET=UC*ET 

EP=UC*EP 
C 
C E-THETA IS CO-POL; E-PHI IS CROSS-POL 
C 


ANG(I)=THETA 
ECV(I)=CABS(ET) 
EX V(I)=CABS(EP) 
ECR=REAL(ET) 
ECI=AIMAG(ET) 
EXR=REAL(EP) 
EXI=AIMAG(EP) 
PHC(I)=ATAN2(ECI,ECR + 1.E-20)/RAD 
PHX(I)=ATAN2(EXI,EXR + 1.E-20)/RAD 
ECX = AMAX1(ECX,ECV(I).EXV(I)) 
500 CONTINUE 
DO 600 I=1,IT 
ECV(I)=AMAXI(ECV(I),1.E-10) 
EXV(I)=AMAXI(EXV(I),1.E-10) 
ECP(I) =(ECV(I)/ECX)**2 
EDP(I) =(EXV(I)/ECX)**2 
ECP(I) =AMAX1(ECP(I)..00001) 
EDP(I)=AMAX1(EDP(I),.00001) 
ECP(I)=10.*ALOG10(ECP(1)) 
EDP(I)=10.*ALOG10(EDP(1)) 
600 CONTINUE 
SIGMA =(ECX**2)*CABS(UC)/(2.*BK) 
SIGDB=10.*ALOG10(SIGMA) 
WRITE(6,*) “SIGMA, IN DB=’,SIGMA,SIGDB 
DO 9000 L=1,IT 
WRITE(8,*) ANG(L),ECV(L) 
9000 CONTINUE 
900 CONTINUE 
STOP 
END 
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201 SUBROUTINE SOLVE(N IES Ue 


202 C 
203 C SUBROUTINE TO SOLVE SYSTEM OF EQUATIONS WITH COMPLEX 
204 COEFFICIENTS. 

205 C CALL ’DECOMP’ FIRST. (FROM MAUTZ AND HARRINGTON) 
206 E 

207 COMPLEX UL(50000),B(500),X(500),SUM 

208 DIMENSION IPS(500) 

209 NPI=N+1 

210 IP=IPS(1) 

211 X(1)=B(IP) 

pile) DO 2 1=2,N 

213 IP =IPS(1) 

214 IPB=IP 

215 IM1=I-1 

216 SUM =(0.,0.) 

217 DO 1 J=1,IM1 

an SUM=SUM + UL(IP)*X(J) 

219 1 IP=IP+N 

220 2 X(I)=B(IPB)-SUM 

221 K2=N*(N-1) 

Np) IP=IPS(N)+K2 

223 X(N)=X(N)/UL(P) 

224 DO 4 IBACK=2,N 

225 I=NP1-IBACK 

226 K2=K2-N 

227 IP] =IPS(I) +K2 

228 IP] =I+1 

229 SUM =(0.,0.) 

230 IP=IPI 

231 DO 3 J=IP1,N 

232 IP=IP+N 

oEe 3 SUM=SUM + UL(IP)*X(J) 

234 4 X(I)=(X(1)-SUM)/UL(IPI) 

235 RETURN 

236 END 

237 SUBROUTINE DECOMP(N,IPS,UL) 

238 C 

239 C SUBROUTINE TO DECOMPOSE SYSTEM OF EQUATIONS. 
240 C FROM MAUTZ AND HARRINGTON. 

241 E 

242 COMPLEX UL(50000),PIVOT,EM 

243 DIMENSION SCL(500),IPS(500) 

244 DO 5 I=1,N 

245 IPS(I)=I 

246 RN=0. 

247 Ji=l 

248 DO 2 J=1,N 

249 ULM =ABS(REAL(UL(J1))) + ABS(AIMAG(UL(J1))) 
250 JI=J1+N 
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C 


IF(RN-ULM) 1,2,2 
1RN=ULM 
2 CONTINUE 
SCL(I)=1./RN 
5 CONTINUE 
NM1=N-1 
K2=0 
DO 17 K=1,NM1 
BIG=0. 
DO 11 I=K,N 
IP=IPS(d) 
IPK=IP+K2 
SIZE=(ABS(REAL(UL(IPK))) + ABS(AIM AG(UL(IPK))))*SCL(UIP) 
IF(SIZE-BIG) 11,11,10 
10 BIG=SIZE 
IPpV =I 
11 CONTINUE 
IF(PV-K) 14,15,14 
14 J=IPS(kK) 
IPS(K)=IPSCUPV) 
IPSUPV)=J 
IS KEP SIPs KZ 
PIVOT =UL(KPP) 
KP1=K+1 
DO 161K FIEN 
KP=KPP 
IP IPS he 
EM =-UL(IP)/PIVOT 
18 ULUP)=-EM 
DO 16 J=KP1,N 
IP=IP+N 
KP=KP+N 
ULdP)= UL(dP) + EM*UL(KP) 
16 CONTINUE 
K2=K2+N 
17 CONTINUE 
RETURN 
END 


SUBROUTINEZMATWW(NWIRES,NW1,NW2,RB,XH, YH,ZH,NT,AT,AT,AK, ZZ) 


C IMPEDANCE ELEMENTS FOR LINEAR BASIS FUNCTIONS. 


S 


COMPLEX CEXP,Z(200),ZZ(15000),CON,CMPLX,EXP 
COMPLEX U0,SUM,SUM1,SUM2,SUM3,SUM4 


DIMENSION ZH(200),XT(128), AT(128).XH(200), YH(200), UU(200) 


DIMENSION D1(200),S1(200),C1(200),ZS1(200) 
DIMENSION XS1(200),YS$1(200) 

DIMENSION CU(200),SU(200) 

INTEGER NT,NWIRES,NW1(4),NW2(4),NS(4) 
PI =3.1415926 
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301 ER=Z--ri 


302 ETA=377. 

303 BK =PI2 

304 U0=(0.,0.) 

305 CON =(0.,1.)*BK*ETA/(4.*P]) 
306 A=AK 

307 € 

308 C DEFINE GEOMETRY TERMS FOR THE WIRE. XH, YH,ZH ARE ALL KNOWNS. 
309 EC 

310 DO 5 L=1,NWIRES 

311 5 NS(L)=NW2(L)-NW1(L)+1 
B12 NS1=NW2(NWIRES)-NW1(1) 
313 NPS=NS1+1 

314 NTRIA=NPS-2 

315 DO 10 N=2,NPS 

316 NO=N-1 

317 I=NW1(1)+N-1 

318 I2=1-1 

319 C 

320 C AVERAGE VALUES 

321 S 

322 ZS1(NO)=.5*(ZH(I) +ZH(12)) 
323 XS1(NO) =.5*(XH(1) + XH(12)) 
324 YS1(NO)=.5*(YH(I) + YH(12)) 
325 DX = XH(1)-XH(12) 

326 DY = YH(I)-YH(12) 

327 D1(NO)=SORT(DX**2+ DY **2) 
328 UU(NO)= ATAN2(DY,DX +1.E-5) 
329 CU(N0)=COS(UU(NO)) 

330 SU(NO) =SIN(UU(NO)) 

331 S1(NO)=DR/D1(NO) 

332 C1(NO)=DZ/D1(NO) 

333 10 CONTINUE 

334 IP=1 

335 WRITE(6,*) *IP=",IP 

336 DO 600 JQ=1,NTRIA 

337 C 

338 C DOING Ii 

339 € 

340 I=IP 

34] J=JO 

342 CC =COS(UU(I)-UU(J)) 

343 TIP=1./D1(I) 

344 TJP=1./Di(J) 

345 T1=D1(1)/2. 

346 T2=D1(J)/2. 

347 SUM =U0 

348 DO 100 K=1,NT 

349 T=T1*XT(K) 

350 Tl=.5+T/D1(I) 
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351 
352 
353 
354 
305 
356 
397 
358 
320 
360 
361 
362 
363 
364 
365 
366 
367 
368 
369 
370 
37 
372 
S873 
374 
375 
376 
S77 
378 
72 
380 
26 | 
BO2 
383 
384 
385 
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387 
388 
389 
520 
39] 
SV 
393 
394 
ENS) 
396 
397 
398 
399 
400 


XI=XS1(1I)+T*CU(D) 
YI=YS1(1I)+T*SU() 
Z1=ZS1(1) 
DO 100 L=1,NT 
TP=T2*XT(L) 
TJ =.5+TP/D1(J) 
XJ=XS1(J)+ TP*CU() 
YJ=YS1(J)+TP*SU(J) 
ZJ=ZS1(J) 
RP =SQRT((XI-XJ)**2 + (YI-YJ)**2 + A**2) 
EXP=CEXP(CMPLX(0.,-RP))/RP 
SUM =SUM + AT(L)*AT(K)*EXP*(TI*TJ*CC-TIP*TJP) 
CONTINUE 
SUM1=SUM*T1*CON*T2 


C DOING [2 


J=JO+1 
CC=COS(UU(I)-UU(J)) 
TIP=1./D1(1) 
Te eA) 
T1=D1(1)/2. 
T2=D1(J)/2. 
SUM =U0 
DO 101 K=1.NT 
T=T1*XT(K) 
TI=.5+T/DI(D) 
XI =XS1(1) +T*CU(1) 
YI=YS1(1)+T*SU(I) 
ZI=ZSI1(1) 
DG 101 L=1,NT 
TP=T2*XT(L) 
TJ =.5-TP/D1(J) 
XJ =XS1(J) + TP*CU(J) 
YJ=YS1(J)+TP*SU(J) 
ZJ=ZS1(J) 
RP=SQRT((XI-XJ)**2 + (YI-YJ)**2 + A**2) 
EXP=CEXP(CMPLX(0.,-RP))/RP 
SUM =SUM + AT(L)*AT(K)*EXP*(TI*TJ*CC-TIP*TJP) 
CONTINUE 
SUM2=SUM*T1*CON*T2 


C DOING J3 


J=IP+1 

J=JQ 
CC=COS(UU(I)-UU(J)) 
TIP=-1./D1(1) 
TJP=1./D1(3) 
T1=D1(1)/2. 


6] 


401 
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T2=D1(J)/2. 
SUM=U0 
DO 102 K=1,NT 
T=T1*XT(K) 
TIl=.5-T/D1(]) 
XI=XS10) + T*CvUud) 
YI=YS1(D+T*SUd) 
ZI=ZSI1 (I) 
DO 102 L=1,NT 
TPZ) 
TJ =.5+TP/D1(J) 
XJ=XS1(J)+TP*Cu() 
YJ=YS1()+TP*SU(J) 
ZJ =ZS1(J) 
RP=SQRT((XI-XJ)**2 +(YI-YJ)**2 + A**2) 
EXP=CEXP(CMPLX(0.,-RP))/RP 
SUM =SUM+AT(L)*AT(K)*EXP*(TI*TJ*CC-TIP*TJP) 
102 CONTINUE 
SUM3=SUM*T1*CON*T2 
C 
C DOING 14 
ce 
J=JQ+1 
CC=COS(UU(I)-UU()) 
TIP=-1./D1(1) 
TJP=-1./D1(J) 
T1=D1(1)/2. 
T2=D1(J)/2. 
SUM=U0 
DO 103 K=1,NT 
T=T1*XT(K) 
TI=.5-T/D1(1) 
AI =XS1(1) + T*CUC) 
Y¥i—Yslijy+ir*sud) 
ZI1=ZS1 (1) 
DO 103 L=1,NT 
TP=T2*XT(L) 
TJ=.5-TP/D1(J) 
XJ = XS1(0J)+TP*Cucy) 
YJ=YS1(J)+TP*SU) 
ZJ=ZS1(J) 
RP=SQRT((XI-XJ)**2 + (Y1-YJ)**2+ A**2) 
EXP=CEXP(CMPLX(0.,-RP))/RP 
SUM=SUM +AT(L)*AT(K)*EXP*(TI*TJ*CC-TIP*TJP) 
103 CONTINUE 
SUM4=SUM*T1*CON*T2 
C 
C IMPEDANCE ELEMENT FOR IP,JQ 
e 
KK=(JQ-1)*NTRIA +IP 
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Z(JQ)=(SUM1+SUM2+SUM3+SUM4) 
600 CONTINUE 
Z(NTRIA)= Z(2) 
Cc 
C FILL THE ENTIRE Z MATRIX USING SYMMETRY PROPERTIES 
C ROW INDEX, I; COL INDEX, J 
Cc 
DO 12 1=1,NTRIA 
DO 12 J=1,NTRIA 
K=(I-1)*NTRIA +J 
ZZ(K)= U0 
WJ =IABS(-J) 
IF(J.GT.NTRIA) GO TO 12 
Ji=U+1 
ZZ(K) = Z(1J 1) 
12 CONTINUE 
CLOSE(2) 
RETURN 
END 
SUBROUTINE PLANEW(NP,XH,YH,ZH,THR,PHR,R) 


PLANE WAVE EXCITATION VECTOR ELEMENTS FOR WIRE AND 
INGIDENGE DIRECTION IS (1 HR.PHR). 
WIRE LIES IN THE X-Y PLANE (Z=0) 


CaP 27a Re 


COMPLEX UO0,C,R(2000),CEXP,EXP,FI1,FI2,S1,DI,CMPLX 
DIMENSION ZH(500),XH(500), YH(500) 
MP2=NP-] 
MT2=NP-2 
U0=(0.,0.) 
CC=COS(THR) 
SS=SIN(THR) 
CP=COS(PHR) 
SP=SIN(PHR) 
UP=SS*CP 
VP=SS*SP 
DO 12 IP=1,MP2 
II=IP 
I=l]+1 
ZS =.5*(ZH(1) + ZH(ID) 
XS =.5*(XH(1) + XH(1)) 
YS=.5*(YH() + YH(ID) 
DX = XH(I)-XH(II) 
DY =YH()-YH(II) 
D1=SQRT(DX**2+DY**2) 
SU=DY/D1 
CU=DX/D1 
C FOR WIRES IN THE XY PLANE SIN(V)=1 AND COS(V)=0 
SV=1.0 
Gy =00 
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C 


C WIRE SEGMENT CALCULATIONS 


G 


c 


60 
61 


62 


A=UP*CU+VP*SU 
B=UP*XS+VP*YS 
C=CMPLX(0.,A) 
EXP=CEXP(CMPLX(0.,B)) 
AA=CC*(CU*CP + SU*SP) 
BB =SU*CP-SP*CU 
PSI=D1*A/2. 
IF(PSI.NE.0.) GO TO 60 
SINC=1. 

GO TO 61 
SINC=SIN(PSI)/PSI 
COSP=COS(PSI) 

FI1 =SINC*D1*EXP/2. 

FI2=(0.,0.) 

IF(ABS(A).LT.1.E-4) GO TO 62 

CSP=COSP-SINC 

IF(ABS(CSP).LT.1.E-4) GO TO 62 

FI2=EXP/C*CSP 
CONTINUE 

SI=FI1 +FI2 

DI=FI1-FI2 


C R-WIRE-THETA 


S 


S 


10 


IF(IP.EQ.MP2) GO TO 10 
RUIP)=AA*SI 
RUP+MT2)=BB*SI 
CONTINUE 


C R-WIRE-PHI 


c 


14 


Irom EO) GOTO 12 
R(UIP-1)=R(P-1) + AA*DI 
R(IP-1 +MT2)=R(P-1 +MT2)+ BB*DI 


12 CONTINUE 


210 RETURN 
END 
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MAb hf hh Hh AA Hh PWWWWWWWWWWNNYNNN YN NY NYY DN RR RR RS Re eS 
Owonaowx7ans WNK OO WAHYAMNHWNK OO WAAIANMAHWNK OW OHA MN Hh WN — © WO CO 


C MAIN PROGRAM *HARLOOP.F* 

C PLANE WAVE SCATTERING FROM A CIRCULAR LOOP IN THE Z PLANE. 

C **** RADAR CROSS SECTION CALCULATION ***** 

C USING HARRINGTON’S FORMULATION FROM THE BOOK ’FIELD COMP. BY 
C MM’ (P.83 TO 95) 


C 


COMPLEX Z(5000),E(250),C(250),EX ,EC,ET,EP,RW(1000),UC, EPHI(500) 
COMPLEX CPHI(500) 
DIMENSION ECP(500),IPS(250), ANG(500),EDP(500),XT(300), AT(300) 
DIMENSION ECV(500),EX V(500),PHC(500),PHX (500) 
DATA PI/3.14159265/ 
DATA IPRINT/1/,ITEST/1/ 
RAD =PI/180. 
ECX=0. 
BK=72-7 Pl 
ETA=377. 
UC=(0,-1)*ETA*BK/(4. *PI) 
PHIRD=0. 
OPEN(1,FILE=’PARAMLST.DAT’) 
READ(1,*) ANGLE,DT,START,STOP,A,B,SEG,AHARR,GCONST,XLOC,XIJPOL 
IPOL=INT(XIPOL) 
LOC=INT(XLOC) 
CLOSE(1) 
NM =INT(AHARR) 
CON=(G677.7BN a2 278K 
OPEN(2,FILE=’OUTGLEG') 
READ(2,*") NT 
DO 1 J=1,NT 
READ(2,*) XT), AT(I) 
CONTINUE 
OPEN(8,FILE=*OUTHARR.DAT'’) 
OPEN(7,FILE= *IHARROUT.DAT’) 
NROW =2*NM +1 
WRITE(6,1300) B,A,NROW NT 


1300 FORMAT(//,5X,’PLANE WAVE SCATTERING BY A CIRCULAR LOOP’,/ 


C 


aD, "USING METHOD IN HARRINGTONS MM TEXT BOOK’,/ 

* xX, LOOP RADIUS (WAVL)=’,F8.4,/,5X,’WIRE RADIUS (WAVL)=’, 

* F8.4,/,5X,, NUMBER OF AZIMUTHAL MODES (INCLUDING ZERO) =",]3, 
* /,5X,, NUMBER OF INTEGRATION POINTS IN PHI=’,14) 


C COMPUTE IMPEDANCE MATRIX ELEMENTS 


C 


C 


WRITE(6,*) "CALLING ZMAT’ 

CALL ZMATWW(NM,A,B,NT,XT,AT,Z) 
WRITE(6,*) "WIRE IMPEDANCE COMPUTED’ 
CALL DECOMP(NROW, IPS, Z) 

WRITE(6,*) °>Z DECOMPOSED’ 


C BEGIN FIELD CALCULATIONS 


C 
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ay 
4 
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74 
75 
76 
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78 
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80 
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82 
83 
84 
85 
86 
87 
88 
89 
90 
a1 
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93 
94 
28) 
96 
oF 
98 
ee, 
100 


PHRO=PHIRD*RAD 
IT=INT((STOP-START)/DT) + 1 
WRITE(7,*) IT,NROW 
DO 500 I=1,]T 
THETA =FLOAT(I-1)*DT+START 
WRITE(6,*) "THETA=', THETA 
THR=THETA*RAD 
PHR=PHRO 
IF(THETA.LT.180.) GO TO 99 
THR =(360.-THETA)*RAD 
PHR=PHRO+PI 
99 CONTINUE 
ET=(0.,0.) 
EP=(0.,0.) 
CALL PLANEW(NM,B,THR,PHR,RW) 
C TRANSMIT VECTOR ELEMENTS ARE TRANSPOSED FORMS OF RECEIVE 
VECTOR 
C ALSO THE THETA COMPONENT GETS A NEGATIVE SIGN 
IF(LOC .EQ. 0) THEN 
IF(IPOL.EQ.1) THEN 


e 
C THETA POLARIZED INCIDENT WAVE (IPOL=1) 
C 
DO 101 L=1,NROW 
101 E(NROW-L-+ 1)=RW(L) 
ESE 
C 
C PHI POLARIZED INCIDENT WAVE (IPOL=2) 
C 
DO 102 L=1,NROW 
102 E(NROW-L+ 1)=RW(L+NROW) 


ENDIF 
WRITE(6,*) ‘CALLING SOLVE’ 
CALL SOLVE(NROW, IPS,Z.E,C) 
WRITE(6,*) "RETURNED FROM SOLVE’ 
DO 210 L=1,NROW 

WRITE(7,*) C(L) 

ET=ET+RW(L)*C(L) 

210 EP=EP+RW(L+NROW)*C(L) 


EESE 
‘3 
CE-THETATS CO*POL E-PHiisS Goss Oe 
C 


DO 221 L=1,NROW 
221  E(NROW-L+1)=RW(L) 
DO 222 L=1,NROW 
oD? EPH (NROW-L+1)=RW(L+NROW)*CEXP(CMPLX(0. ,PI/2.)) 
CALL SOLVE(NROW, IPS, Z,E,C) 
CALL SOLVE(NROW, IPS, Z,EPHI,CPHI) 
DO 220 L=1,NROW 
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147 
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220 


WRITE(7 :*) C(_) + CPHKL) 
ET=ET+RW(L)*C(L)+RW(L)*CPHI(L) 
EP=EP+RW(L+NROW)*C(L)+RW(L+NROW)*CPHI(L) 

ENDIF 

EC—UG*ET 

ExX=UG-EP 

ANG(I)=THETA 

ECV(I)=CABS(EC) 

EXV(I)=CABS(EX) 

ECR=REAL(EC) 

ECI = AIMAG(EC) 

EXR=REAL(EX) 

EXI=AIMAG(EX) 

PHC(1)=ATAN2(ECI,ECR + 1.E-20)/RAD 

PHX(I)=ATAN2(EXI,EXR+1.E-20)/RAD 

ECX = AMAX1(ECX,ECV(]),EXV(I)) 


500 CONTINUE 


DO 600 1=1,IT 
ECV(I)=AMAX1(ECV(1),1.E-10) 
EXV(I)=AMAXI(EXV(I),1.E-10) 
ECP(I)=(ECV(I)/ECX)**2 
EDP(I) =(EXV(I)/ECX)**2 
ECP(I)= AMAX1(ECP(I),.00001) 
EDP(I)=AMAX1(EDP(I),.00001) 
ECP(I) = 10.*ALOG10(ECP(I)) 
EDP(I)=10.*ALOG10(EDP(1)) 


600 CONTINUE 


c 


SIGMA =(ECX**2)*4.*P] 
SIGDB=10.*ALOG10(SIGMA) 

WRITE(6,*) "BACKSCATTER, IN DB=’,SIGMA,SIGDB 
OPEN(2,FILE=’RCS.DAT’) 

WRITE(2,*) SIGMA,SIGDB 

CLOSE(2) 


C PRINIOPIELD POINTS 


c 


DO 9000 L=1,IT 
WRITE(8,*) ANG(L),ECV(L) 


5016 FORMAT(5X,F6.1,3X,2(F8.4,3X,F7.1,3X,F7.2,3X)) 
9000 CONTINUE 

900 CONTINUE 

9998 STOP 


END 

SUBROUTINE SOLVE(N,IPS,UL,B,X) 
COMPLEX UL(5000),B(250),X(250),SUM 
DIMENSION IPS(250) 

NP1=N+1 

IP=IPS(1) 

X(1)=B(IP) 

DO 2 I=2,N 
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IP=IPS(I) 


IPB=IP 
IM1=I-1 

SUM =(0.,0.) 
DO 1 J=1,IM1 


SUM =SUM + UL(P)*X(J) 


1 IP=IP+N 
2 X(1)=B(IPB)-SUM 


K2=N*(N-1) 
IP=IPS(N) +K2 
X(N)=X(N)/UL(IP) 
DO 4 IBACK=2,N 
I=NP1-IBACK 
K2=K2-N 
IPI=IPS()+K2 
IPl=I+1 

SUM =(0.,0.) 
hel 

DO 37 —l1F iN 
IP=IP+N 


3 SUM=SUM + UL(IP)*X(J) 
4 X(1) =(X(1)-SUM)/UL@PI) 


RETURN 

END 

SUBROUTINE DECOMP(N,IPS,UL) 
COMPLEX UL(5000),PIVOT,EM 
DIMENSION SCL(250),1PS(250) 
DO5 [=1.N 

IFS()=! 

RN =0. 

ipl 

DO =I1N 

ULM =ABS(REAL(UL(J1))) + ABS(AIMAG(UL(1))) 
JI=J1+N 

TRGEN- WIEN hse 


1] RN=ULM 
2 CONTINUE 


SCL(I)=1./RN 


5S CONTINUE 


NM1=N-1 

K2=0 

DO 17 K=1,NM1 

BIG=0. 

DO 11 I=K,N 

IP=IPS(1) 

IPK=IP+K2 
SIZE=(ABS(REAL(UL(IPK))) + ABS(AIMAG(UL(IPK))))*SCL(IP) 
IF(SIZE-BIG) 11,11,10 


10 BIG=SIZE 


IPV=] 
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11 CONTINUE 
IF(IPV-K) 14,15,14 
14 J=IPS(K) 
IPS(K)=IPSUPV) 
IPSUPV)=J 
15 KPP=IPS(K)+K2 
PIVOT =UL(KPP) 
KP1=K+1 
DO 16 I=KP1,N 
KP=KPP 
IP=IPS(1)+K2 
EM =-UL(IP)/PIVOT 
18 UL(IP)=-EM 
DO 16 J=KPI1,N 
IP=IP+N 
KP=KP+N 
ULUP)=ULGUP)+EM*UL(KP) 
16 CONTINUE 
K2=K2+N 
17 CONTINUE 
RETURN 
END 
SUBROUTINE ZMATWW(NM,A,B,NT,AT,AT,Z) 
Cc 
C *** MODS FOR LOOP -- USING HARRINGTON’S TEXT BOOK EQUATIONS AS A 
C CHECK FOR OTHER METHODS. NM IS THE NUMBER OF AZIMUTHAL MODES 
USED 
Cc 
COMPLEX Z(5000),CON,CMPLX,FK 
DIMENSION XT(300),AT(300) 
Pi 391415926 
BR=2- Fl 
CON=CMPLX(0. ,PI*377.*BK*B) 
NROW =2.*NM +1 
DO 10 I=-NM,NM 
DO 10 J=-NM,NM 
J=J+NM+1+(1+NM)*NROW 
Z(1J)=(0.,0.) 
10 CONTINUE 
C 
C ONLY DIAGONAL ELEMENTS ARE NONZERO. ALTHOUGH SYMMETRY EXISTS 
BETWEEN 
C Z(-N,-N) AND Z(N,N) IT IS NOT BEING USED. 
Cc 
DO 20 I=-NM,NM 
J=1+NM+1+4+(1+NM)*NROW 
JIP=]+1]1 
IM=I-1 
Z(J)=(.5*FKUM,B,A,NT,XT,AT)+.5*FK(IP,B,A,NT,XT,AT)- 
7 “A Bie Ey=<2*F KB ANT, X EP, AL))*CON 
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20 CONTINUE 
RETURN 
END 
SUBROUTINE PLANEW(N,B,THR,PHR,R) 


PLANE WAVE RECEIVE VECTOR ELEMENTS FOR WIRE USING THE 
FORMULATION FROM HARRINGTON’S BOOK. N IS THE NUMBER OF 
AZIMUTHAL MODES. NOTE THAT B(N)=R(-N) (B IS EXCITATION AND 
R IS RECEIVE). 


OOO @Oy@ 


COMPLEX R(1000), CEXP,EXP,CMPLX 
PI=3.14159 
BK=2.*PI 
CT=COS(THR) 
ST =SIN(THR) 
RR=BK*B*ST 
NROW =2*N+1 
C DO THETA RECEIVE COMPONENTS FIRST 
DO 10 I=-N,N 
IP=I+1 
IM =I-] 
EXP=CEXP(CMPEX(O31-PHR)) 
RUI+N+1)=-PI*B*(0.,1.)**I*EXP*(BESSJ(IP,RR)+BESSJ(IM,RR))*CT 
10 CONTINUE 
C NOW DO PHI RECEIVE COMPONENTS 
DO 20 I=-N,N 
IP=1+1 
IM =I-1 
EXP=CEXP(CMPLX(0.,I*PHR)) 
R(I+N+1+NROW)=PI*B*(0.,1.)**(1+ 1)*EXP*(BESSJ(IP,RR) 
* -BESSJ(IM,RR)) 
20 CONTINUE 
RETURN 
END 
COMPLEX FUNCTION FK(N,B,A,NT,XT,AT) 
COMPLEX SUM,EXP1,EXP2,CEXP,CMPLX 
DIMENSION XT(300),AT(300) 
PI=3.14159 
BK =2.*PI 
CC=1./BK 
P1=(2.*PI-0.)/2. 
P2=—(2. Pl +0 )2. 
SUM =(0.,0.) 
DO 10 1=1,NT 
PHI=P1*XT(I)+P2 
RR=SQRT(4. *SIN(PHI/2.)**2 + (A/B)**2) 
EXP1=CEXP(CMPLX(0.,-BK*B*RR)) 
EXP2=CEXP(CMPLX(0.,-N*PHI)) 
SUM=SUM+AT(I)*EXP1*EXP2/RR 
10 CONTINUE 
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FK=SUM*P1*CC 
RETURN 

END 

FUNCTION BESSJ(NN,X) 


C KETURNS THE BESSEL FUNCIION BOF ORDER N (> 1) AND REAL 
C ARGUMENT xX. 


PARAMETER (IACC=40,BIGNO=1.E10,BIGNI = 1.E-10) 
IF(NN.LT.0) N=-NN 
IF(NN.GE.0) N=NN 

KC=3 

IF(N.EQ.0) KC=1 

IF(N.EQ.1) KC=2 

GO TO (1,2,3),KC 

BESS) =BESSJO(X) 

GOTO 4 

BESS) = BESSJ1(X) 

Come. 

BESSJ=0. 
IF(ABS(X).LT.1.E-5) GO TO 4 
TOX =e 
IF(X.GT.FLOAT(N)) THEN 

BJM =BESSJ0(X) 

BJ =BESSJ1(X) 

DO 11 J=1,N-1 
BJP=J*TOX*BJ-BJM 
BJM=BJ 
BJ=BJP 

11 CONTINUE 

BESSJ=BJ 

ELSE 

M =2*((N+INT(SQRT(FLOAT(IACC*#N))))/2) 

BESSJ=0. 

JSUM =0. 

SUM =0. 

BJP=0. 

BJ=1. 

DO 12 J=M,1,-1 
BJM =J*TOX*BJ-BJP 
BJP=BJ 
BJ=BJM 
IF(ABS(BJ).GT.BIGNO) THEN 

BJ =BJ*BIGNI 
BJP=BJP*BIGNI 
BESS] =BESSJ*BIGNI 
SUM =SUM4BIGNI 
ENDIF 
IF(JSUM.NE.0) SUM=SUM +B] 
JSUM=1-JSUM 
IF(J.EQ.N) BESSJ=BJP 
12 CONTINUE 
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376 
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380 
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SUM =2.*SUM-BJ 
BESSJ=BESSJ/SUM 

ENDIF 

CONTINUE 

IF(NN.LT.0) BESSJ =(-1.)**N*BESS]J 

RETURN 

END 

FUNCTION BESSJO0(X) 


C BESSEL FUNCTION OF 0 ORDER, REAL ARGUMENT X 
C (SEE “NUMERICAMPRECIPES? Ply 2) 


Cc 


C 


REAL*8 Y,P1,P2,P3,P4,P5,Q1,Q2,Q3,Q4,Q5,R1,R2,R3,R4,R5,R6, 
* $1,82,83,84,S5,86 

DATA P1,P2,P3,P4,P5/1.D0,-.109862827D-2,.2734510407D-4, 

* -.2073370639D-5,.2093887211D-6/ 

DATA Q1,Q2,Q3,04,Q5/-.1562499995D-1,.1430488765D-3, 

* -.6911147651D-5,.7621095161D-6,-.934945152D-7/ 

DATA R1,R2,R3,R4,R5,R6/57568490574.D0,-13362590354.D0, 

* 651619640.7D0,-11214424.18D0,77392.33017D0,-184.9052456D0/ 
DATA S1,S2,$3,84,85,86/57568490411.D0, 1029532985.D0, 

* 9494680.718D0,59272.64853D0,267.8532712D0, 1.D0/ 
IF(ABS(X).LT.8.) THEN 

Y=X**2 

BESSJO=(R1 + Y*(R2+ Y*(R3+ Y*(R4+ Y*(R5 + Y*R6)))))/ 

* (SIF Y*(S2+Y*(S3+Y*(S4+ Y*(S5 + Y*S6))))) 

ELSE 

AX =ABS(X) 

Z=8./AX 

Y=Z**2 

XX = AX-.785398164 
BESSJO=SORT(.636619772/AX)*(COS(XX)*(P1 + Y*(P2 + Y*(P3 + 
* Y*(P4 + Y*P5))))-Z*SIN(XX)*(QI + Y*(Q2+ Y*(Q3 + 
 OYS(Os ¥*O5))))) 

ENDIF 

RETURN 

END 

FUNCTION BESSJ1(X) 


C BESSEL FUNCTION B OF ORDER 1, REAL ARGUMENT X 
C GEE “NUMERICAL RECIPES’, P.173} 


c 


REAL*8 Y,P1,P2,P3,P4,P5,Q1,Q2,Q3,Q4,Q5,R1,R2,R3,R4,R5,R6, 
* §1,S2,S3,84,S5,S6 

DATA P1,P2,P3,P4,P5/1.D0,.183105D-2,-.3516396496D-4, 

* 2457520174D-5,-.20337019D-6/ 

DATA Q1,Q2,Q3,04,Q5/.04687499995D0,-.2002690873D-3, 

* 8449199096D-5,-.99228987D-6,.105787412D-6/ 

DATA R1,R2,R3,R4,R5,R6/723626 14232. D0,-7895059235.D0, 

* 242396853. 1D0,-297261 1.439D0, 15704.4826D0,-30. 16036606D0/ 


i 


401 
402 
403 
404 
405 
406 
407 
408 
409 
410 
411 
412 
413 
414 
415 
416 
417 


DATA S1,S82,S83,S84,S5,86/144725228442.D0,2300535178.D0, 
* 18583304.74D0,99447.43394D0,376.999 1397D0,1.D0/ 
IF(ABS(X).LT.8.) THEN 

Y=X**2 

BESSJ1=X*(R1 + Y*(R2+ Y*(R3 + Y*(R4+ Y*(RS5 + Y *R6)))))/ 
* (S1+Y*(S2+ Y*(S3 + Y*(S4+ Y*(S5 + Y*S6))))) 

BESE 

AX =ABS(X) 

Z=8./AX 

Y= 

XX =AX-2.356194491 
BESSJ1=SQRT(.636619772/AX)*(COS(XX)*(P1+ Y*(P2+ Y*(P3+ 
* Y*(P4+ Y*PS5))))-Z*SIN(XX)*(Q1 + Y*(Q24+ Y*(Q3 + 

* Y*(Q44+ Y*Q5)))))*SIGN(1.,X) 

ENDIF 

RETURN 

END 
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ADDITIONAL PLOTS OF CURRENT AND RMS ERROR 


Error versus Segment Size 


Circular Polanzation 
0.4 - Loop Radius =0.1 


ae | Wire Radius=0.001 
oa Ng lambda = 320 ° 


0.34 Incidence Angle = 20deg 


0 0.02 004 006 0.08 0.1 0.12 O14 O16 O18 0.2 
Segment Size. Wavelengths 


Error versus Segment Size 





| Circular Polanzation 
0.4 + Loop Radtus=0.1 


Wire Radius=0.001 
| Ng lambda = 320 
03= Incidence Angle = 40deg ° 


0 0.02 004 006 0.08 0.1 012 014 O16 #018 0.2 
Segment Size, Wavelengths 


7 


A 
i 


Normahzed RMS Error 


Normahzed RMS ftrror 





Error versus Segment Size 







Circular Polanzation 
Loop Radtus=0. 1 

Wire Radius=0.001 
Ne/lambda= 320 
Incidence Angle =60deg 





0.02 0.04 0.06 0.08 0.1 Oere 0.14 0.16 0.18 0.2 


Segment Size, Wavelengths 


Error versus Segment Size 


Circular Polanzation | 
Loop Radtus=0.1 


Wire Radius=0.001 
Ng lambda = 320 | 
Incidence Angle =80deg 





0.02 0.04 0.06 0.08 0.1 0.12 0.14 0.16 0.18 0.2 


Segment Size, Wavelengths 


75 


Normalized RMS Error 


Normalized RMS Error 


Error versus Segment Size 





Ph; Polarization 
Loop Radius=0.1 
Wire Radius=0.001 
Ng/lambda=320 
Incidence Angle =Odeg e 
° 
ie) 
1@] . + 
g one a + 


t 


dus. O04 006 O08. 0.1 SminO™ Guin cuG) mnROUiQn EET 


Segment Size, Wavelengths 


Error versus Segment Size 


Phi Polanzation 

Loop Radius=0.1 

Wire Radius =0.001 

Ng lambda = 320 
Incidence Angle =20deg 





0.02 0.04 0.06 0.08 0.1 Os). 0.14 0.16 0.18 0.2 


Segment Size, Wavelengths 


76 


Normalhzed RMS Error 


Normalized RMS Error 


Error versus Segment Size 





| Phi Polanzation 
0.4 Loop Radtus=0. 1 


Wire Radius=0.001 
| Ng lambda = 320 
0.3 | Incidence Angle =40deg 
| 


0.08 - = . 4 


ys eo Da a rc 
° 0.02 0.04 0.06 0.08 0.1 CS Be 0.14 0.16 0.18 0.2 


Segment Size. Wavelengths 


Error versus Segment Size 
0.5 ae oe een ee 


Phi Polarization 
0.4 - Loop Radius=0.1 


Wire Radius=0.001 
Ne lambda = 320 
‘ae Incidence Angle =60deg 


© 
4. 
rn 
I 
j ern) 


eee eee 


om) 
> ue 
— tn 


© 
© 
in 





© 


0.02 0.04 0.06 0.08 0.1 Onhe 0.14 0.16 0.18 0.2 


Segment Size, Wavelengths 


1a) 


Current, amperes 


Phase, degrees 


x107 Current Magnitude vs Loop angle 





\ / Circular Pol. 
i 4 ; Segment size =0.15 
P F Loop Radius=0. 1 
s we: Wire Radius=0.001 
0.5 Pal Wire B 
4 ome Inc. Angle =60deg | 
0 50 100 150 200 0 300 350 
Loop Angle. Phi 
Phase vs Loop angle 
150 = | 
| 
| ‘ | 
100; 4 
") | 
| 
50 - | 4 
<< | 


eee | euclat Pol. ae 
100+ . gment size =0. 


Loop Radius=0.1 





| Wire Radius=0.001 
150+ Nga 4s 

) Inc. Angle =60deg 

0 50 100 150 200 250 300 350 


Loop Angle, Phi 


78 


Current, amperes 


Phase, degrees 


x07 Current Magnitude vs Loop angle 





Cron 
nl 
3,86 . ; ; 
| Ue oO °o 
3+ ; ee i 
| /- 
UeTe / 
pas 
ug P y 
15+ a 
| Circular Pol. | 
| 3 Segment size =0.1 | 
foe 
. Loop Radius=0.1 | 
, Wire Radius=0.001 
ee Ne=32 = 
; Inc. Angle =40deg | 
a rr St 8 Ee ee 
0 50 100 150 200 25) 300 350 
Loop Angle, Phi 
Fhase vs Loop angle 
S00 = 2 
150 - 2 | 
| 
100 - 
50 - 4 


eye bse 


-100 —— | Circular Pol. 


oc Segment size =0.1 
Loop Radius=0.1 





oe NN Wire Radius =0.001 
| Ne=32 

-200 - Inc. Angle =40deg | 
0 30 100 150 200 250 300 350 


Loop Angle, Phi 


is 


Current, amperes 


Phase, degrees 


x10- Current Magnitude vs Loop angle 

















4.5 - 
ie 
‘ 5 I po 
ie 
25 : 
vA 2 
ae / 
| “ 
is : ve 
i Circular Pol. 
Segment size =0.2 | 
oe Fe Loop Radius =0.1 ) 
0.5 - ee - 
F Inc. Angle =40deg | 
0 30 100 150 200 250 300 350 
Loop Angle. Phi 
Phase vs Loop angle 
aK = 
150 - 2) 
100'= ; 
30 - z 
Os _ 
-SQ - 4 
3 | Circular Pol. 
me <4 | Segment size =0.2 
| | Loop Radius=0. 1 
2504 | 
eg \ Wire Radius=0.001 
4  Ng=64 
-200 - Inc. Angle =40deg 
0 30 100 150 200 250 300 350 
Loop Angle. Phi 


80 


Current, amperes 


Phase, degrees 


Current Magnitude vs Loop angle 


x10 ° 










FOS 
\- 
Phi Pol. / 
Te Ys , _ Segmenm size =0.2 
. es Loop Radius =0.2 
4 Kise Radius =0.001 
‘ Inc. Angle =40deg 
0 30 100 150 -00 2a) 300 350 
Loop Angle. Phi 
Phase vs Loop angle 
| 
o tk. eee 
Ly = = 
| 
$0) - : 4 
U- ° SO sco te ° | 
50 - 8 oO 
| Phi Pol. 
-100 + Segment size =0.2 


Loop Radius=0.2 | 
Wire Radius=0.001 | 


Sei Ng=64 
Inc. Angle =40deg | 


0 50 100 130 200 250) 300 350 
Loop Angle, Phi 


8] 


Current, amperes 


Phase, degrees 


Current Magnitude vs Loop angle 





; \ E| ane | + eee 


i \ 
i \ | | larg Phi Pol | 
: \ | | | Segment size =0. / 
E a | } Laop Radius=0. 
Os> AE a \} te ie | * 
7 / / Wire Radius=0. 
| 2 : Ng=32 
; Inc. Angle = 40de g 
0 50 100 150 200 250 300 350 
Loop Angle. Phi 
Phase vs Loop angle 
| 
| 
oe ee we ——— : 
100 a 3 8 3 5 oe Fo, . 5 85 8 ® 4 
. / | o | 
| | 
50 - | } 1 
- | : 
3 : b Pra)! 


Z | i | 
-50 - 6 y \ 
) as ¢ 555 igi Pol. 
gment size =0.1 


Loop Radius=0.5 | 





-100 - Wire Radius =0.001 
: Ng=32 
Inc. Angle =40deg 
0 50 100 150 200 256 300 350 
Loop Angle. Phi 


82 


INITIAL DISTRIBUTION LIST 


Library, Code 52 
Naval Postgraduate School 
Monterey, CA 93943-5002 


Defense Technical Information Center 
Cameron Station 
Alexandria, VA 22304-6145 


Chairman, Code EC 

Department of Electrical and Computer Engineering 
Naval Postgraduate School 

Monterey, CA 93943-5002 


Commander 

Naval Air Test Center 

SY90 Mr. L. Krouse 

Patuxent River, MD 20670-5304 


Commander 

Naval Air Test Center 

SY95BW Mr. B. Walter 

Patuxent River, MD 20670-5304 


Professor D. C. Jenn, Code EC/Jn 

Department of Electrical and Computer Engineering 
Naval Postgraduate School 

Monterey, CA 93943-5002 


Michael A. Morgan, Code EC/Mw 

Department of Electrical and Computer Engineering 
Naval Postgraduate School 

Monterey, CA 93943-5002 


83 


Gye 1S 














Thesis 
W22415 
c.l 


Walter 


The use of conformal 
subdomain basis functions 
in the method of moments 


computations for a thin 
wire. 





77 


= = & seu “ee 


pietsut. 12 
+ oom, oe 


"yar ae oe 
=—" 
ow 
epee Ban w: 
ae 
qm eur ate we 
teen mw 


eo 
eon = « 
Se a eae aa one? 

te a ao « 
P - » =o” 
R Cy wae? 

. * 


a 
neo et ormteves™ o = 
r , 


aa 
ar « 


i 
ace t® 


- At onl? 
im * 


il 


| 
| 


= @5% 
omnes 

gear ene ewee 
” . . 


\\ 


| 
| 
07 9 


pucane te Y a 
= quan tee fa 
gn ea dee OI = ge 
ey 


| 
| 
| 


| 
| 


i 


i 


| 


- 4 
q tee oar 
yoke ow ee 


| 


| 


pepe Dad stat . : 
»~ Wee ¥, pugen woe ETRE” 4 ,. o sé Pe 
ue o ¥ * Z 
- .e* "t ae eo 8 
o eytatsty 
e See 


| 


| 


ogee nooo ate 


an regres te’ owen 


pare ee 
owe 


pea me" Fe * 
. . ~oeare’ 
' = aiert = 
~ enn wane aes Sawree 


| 


| 
| 
3 2768 000359 


| 


re nated) 
exe yow Os & ae 
io re ae Oe 


| 
| 


| 


-7e 2 
«meee 


| 
| 


ac 
< 
jan 
2 
—l 
x 
O= ——— 
Te = 
Y 
if 
— 
QO 
=) 
QO 


Hil 
ll 


| 


rere ft vee out y 
agen ee . ree 
: Papen rr eey Fa) coax! : 
care yeiets sere jieree’ 
yeyan a2 on" 


oa 


| 
| 


| 


pig yenegee eee 

hy verre - ye 

> ee . : ame" - 4 

nat To. a , “ ; o F; “se 

on ete : ; se fee 
Vnae = L ae S 5 E canes ere 

ona : ba o eta oP cp aan saree? 

7 ae ‘ : woh bs ras 009 EP © 

pT mie” 

ye BE at cas coe? 


rane come” 
we aperahege nee 
one yet ge aoe en fe 
rarene’ see 1 aperaam me carat ores? 
i: i nae Meee act Nj. conen owen oe 
geese CS ie eta pore er 


. We 0 bev e orane se 
coma Fe bd = . - ETP Far enleae 
=o agen 

7 - eae ape 


- rp. o 

erart 

* 
8a 
Fi . Pes hnend com Tato 
ied #: Go 

hte = OR ree ee 
age zerswane® ane = a Se pac euns meme 
ee ath POT er nis sap oe tee 


wa eneene Sanh 
‘gegugis bPrrs* b 
cree eT TT alain a™ K 
ag gees 08 yute®=? ogee ‘ 4 : 
a ~. « 
"4 “ rep nner Sone ’ sonore 
sell GPE ED b x pie faatetetgemmitsn ae 
Re Ph > = no eer oe 
net To ene te 
¥ pa ee eet. 


") s 
prea e Prem pre eR 
ean ea tePe ane gba ane > 
<a apts =" preg co eee” ayre® 
Tqape ore a? Catal pepe 
etn baie Pr 


res 
gota ® 
nee 


Pts 
Saat te ©! 
>. 


Siete 


sae eeaetnece 


wee 





ets 
om 
eee 
iat dr omen 
rye S 
Sagi tegtinar8 Ff ot ood ei cgae> 
et 
Pevns 
siatigenay 
> 8 cas : Payee 
~ a @ Th0%. : rede Tyeae eer adh 
, : Met prha cer dane 
5 atari seen te 
~ eniar ee OT kde 
= o* pert ra 
pa 


